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ABSTRACT 

DESAP 2 is a finite element program for computer-automated, minimum weight design of 
elastic structures with constraints on stresses (including local instability criteria) and buckling 
loads. No limits are placed on the number of load conditions for stress-constrained design, but 
only one of these load conditions can be chosen as the potential buckling load. A substantial nortior] 
of DESAP 2, particularly the analysis of the prebuckling state, is derived from the SOLID SAP 
finite element program developed at the University of California, Berkeley. The stress -con r- 
strained design is based on the classical stress ratio method, which drives the design towards a 
fully stressed state. The constraints on the buckling load are handled by solving the appropriate 
optimality criterion by successive iterations. During each iteration, the element sizes deter - 
mined by the stress ratio method are used as the minimum size constraints. The element subrou-] 
tines have been organized in a manner that permits the user to make additions and changes with a 
minimal programming effort. Consequently, DESAP 2 can readily be changed into a special- 
purpose program to handle the user’s specific design requirements and failure criteria. 

DESAP 2 is a companion program of DESAP 1: ”A Structural Design Program with Stress 
and Displacement Constraints. " With the exception of a few cards, the same input data deck 
can be used for both programs. 

This is Volume I of three volumes. 
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A.l 


A. INTRODUCTION 


DESAP 2 is a finite element program for automated design 
(synthesis) of linear-elastic structures under static loads. The 
design objective is to find the element sizes (cross-sectional areas, 
plate thicknesses, etc.) that minimize the total structural weight. 

The layout of the structure is not changed during the design procedure. 

primary constraints used in the synthesis algorithm are 
upper limits on stresses and lower bounds on buckling loads. The 
stress limits may be prescribed in the form of yield criteria, local 
instability criteria, or both. There are no restrictions on the 
number of load conditions that may be imposed on the structure for 
the stress-constrained design. Considerations of economy, however, 
limit the number of load conditions for buckling-constrained design 
to one. 

The program also allows the use of secondary constraints , which 
consist of minimum allowable element sizes and size proportion con- 
straints (the requirement that the sizes of specified elements be 
equal, or have a prescribed ratio). 

The method of design is very similar to that used for stress 
and displacement constraints in DESAP 1[1]; in fact, many parts of 
DESAP 1 and DESAP 2 are identical. The design procedure is an 
iterative process, each iteration consisting of four parts: 

1) Prebuckling analysis of the current design. 

2) Redesign with respect to stress constraints based on the 
results of prebuckling analysis. 



3) Buckling analysis of the current design under the action of 
the internal forces obtained in part 1). 

4) Redesign with respect to buckling; the sizes of elements de- 
termined in part 2) are used as the minimum size constraints. 

The entire computer program is built around the SOLID SAP finite 
element program developed by E. L. Wilson[2], It was necessary, of 
course, to carry out extensive modifications of the existing SOLID SAP 
subroutines in order to accomodate the special requirements of the 
redesign operations. Apart from the modifications, several major addi- 
tions to the program were made. These included the entire buckling 
analysis package, all the redesign subroutines, and the subroutines for 
the shear panel element. 

The classical stress ratio method is employed in the redesign with 
respect to stress constraints . This procedure will drive the final 
design to the fully stressed design, which does not necessarily coincide 
with the minimum weight distribution of the material. Unfortunately, 
at the present state-of-the-art, more rigorous design methods are en- 
tirely impractical for structures of reasonable size, since they re- 
quire the inversion of the structural stiffness matrix, or the use of 
numerical search techniques. 

For the constraints on the general buckling loads , the redesign 
procedure described in [3,4] is used. The redesign formulas are 
derived directly from the optimality criterion; consequently, the 
element sizes are driven towards the minimum weight design. 

The design procedure for stress as well as the buckling constraints 
is based on the assumption that the loading is independent of the 
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element sizes (dead loading of prescribed magnitude). Allowance has 
been made for size-dependent loads,, e.g., thermal stresses and gravity 
loading, but these loads must be used with special precautions. 

As pointed out repeatedly in existing literature, neither the fully 
stressed design, nor the optimal weight design are necessarily unique, 
i.e., they have a local rather than a global character. It follows, 
therefore, that the choice of the initial design plays an important role 
in determining the design to which the synthesis algorithm converges. 
Numerous examples seem to indicate, however, that the weight differences 
between the various converged designs are slight, although there may be 
large differences in the distribution of the material. 

A very important feature of the program is the organization of 
the element subroutines such that they can easily be adapted to the 
user's special requirements with the smallest possible programming 
changes. In making this provision, we recognized the fact that it is 
virtually impossible to create an all-purpose synthesis program. The 
reason is that structural design requires a much larger and more varied 
volume of input information than analysis. In the case of beam elements, 
for example, the cross-sectional area, the two moments of inertia and 
the torsional constant suffice for the purposes of analysis. In de- 
sign we must add the section moduli, local buckling information and the 
yield criterion; in addition we must specify how all these properties 
vary with the design variable. Since each design situation may use 
elements of special construction and shape, and different design cri- 
teria, it is clearly impractical to make provisions for all the 
possible forms of design data in a single progreim. 


IL 


Finally, it must be pointed out that DESAP 2, despite its flexi- 
bility, can still handle only a limited amount of design information. 
Consequently, the results of the program are to be taken as a pre- 
liminary design in the sense that it gives the desired proportions of 
the structure, but detailing must be still carried out by the designer 
through conventional design practices. 
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B, THEORETICAL BACKGROUND 
B. 1 Element Properties 

The size of a typical (i^^) element is denoted by A^. It repre- 
sents the cross-sectional area or the panel thickness for one and two- 
dimensional elements, respectively. The weight of the element can be 
written as 


Wi = p^A. , i = 1, 2, I (B.1.1) 

where is the unit weight of the element. If the minimum cost, rather 
than weight is the design objective, then p^ should be taken as the 
unit cost. 

It is seldom desirable to have each as an independent 
design variable. Equal size constraints can be imposed by introducing 
the design variables D^, m = 1, 2, . . . , M, where M £ I, which are in- 
dependently variable, and expressing each element size in the form 

A. = n.D , (B.1.2) 

where is called the design variable fraction of the element. Both 
and m must be specified for each element. Equal sizing of two 
elements is obtained, for example, by prescribing the same m for each 
element together with = 1- In addition^ the scheme permits propor- 
tional sizing if the same m, but different values of x)^ are used. 

The stiffness matrix of each element is restricted to the form 

[K.] =[kf^5]A. + [kf2)]A"i 


(B.1.3) 
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where [K^] is the element stiffness matrix, [kP^] and [kP^] are the 
unit stiffness matrices of the element, and is the inertia exponent . 
The first part of (B.1.3) represents the action of direct stresses, 
whereas the second part is due to bending or torsion. The value of 
depends on the physical nature of the design variable. For example, 
n^ = 1 for thin-walled beams if the wall thickness only is being varied. 
If all dimensions of the cross section are scaled uniformly, 
then n^ = 2. For plates, where the thickness is subject to design, 

we have n^ = 3 . The unit stiffness matrices of each element are stored 
separately on an external storage device, together with m, and n^. 

It should be noted that n^ is determined by the relationship between 
the gize and the moment of the inertia of the element: 


I. 

1 


n. 

j.A.^ 

•^1 1 


(B.1.4) 


where j. is the unit moment of inertia. 

■> 2 . 

The vector of internal forces {N^} of an element is recovered from the 
element nodal displacements {u^} by the formula 


{N.} = [S.]{u.} + {T.} 


(B.1.5) 


If the element stiffness matrix is given by (B.1.3), then the force 
recovery matrix [Sj^] has the same form: 


[S.] = [s{^^]A. + . 


(B.1.6) 


where [s^^^] and [sP^] are the unit force recovery matrices . The 
force vector {T^} also consists of two parts: 

{T.} = {tf^^A, + {tf°h . 

1 111 


(B.1.7) 
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{t } being the unit thermal force vector and {t } represents the contri- 
bution of size-independent loading; e.g., fixed-end forces of beam ele- 
ments due to dead loading. Again, [sP^], [sP^], {t^^^} and {t } 
are stored for each element. 

The element load vector {Q^^} is similarly separated into two components: 
(Qi> = , (B.1.8) 

where is the unit load vector due to size-dependent loading (thermal 

and gravity loads), and represents the contribution of dead loads; 

both vectors are stored for each element. 

The form of equations (B,l,7) and (B.1.8) allows only for a uniform 
temperature increase within an element, i.e.,, it as.sumes that thermal 
expansion causes extension without bending. This excludes, for example, 
the effects of thermal gradients through the thickness of plates, for which we 
would require additional terms of the type {t^ ^ 
for {T^} and {Q^}, respectively. 

Up to now the discussion has been confined to the properties of the 
elements associated with linear analysis, i.e. the analysis of the pre- 
buckling state. In the analysis of general buckling, the elastic stiff- 
ness matrix [K^] of each element must be supplemented by its geometric 
stiffness matrix [G^] . The latter represents the contribution of the 
nonlinear terms in the strain-displacement relations to the strain energy 
of the element, the linear terms being accounted for by [K^] . 
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The buckling analysis of DESAP 2 is limited to bifurcation in- 
stability (classical, Euler-type buckling), in which case only the 
quadratic terms in the components of the rotation vector {00} are in- 
cluded in the formation of [G^] • The contribution of these terms to 
the strain energy of an element is see Ref. [5], p. 7.1-6. 


= T-^[Ca °-Ha 0)0) 2 . (a O.a 0)0) 2 ^ ^ 0^ 0^ 2 


(B.1.9) 


- 0 .0 
-2t 0)0) - 2t 0)0) 

xy X y yz y z 


2 T 0) 0) ]dV, 
zx z X 


where the superscript "0" denotes that the stresses are evaluated in 
the prebuckling state, and V refers to the volume of the element. 

The rotations are given in terms of the buckling displacements u,v,w 


by 


1 ,9w Sv. 1 .3u 3w., 

2^3y “ 3z^’ ‘^^y " 2^3z 3x^ ’ 


1 _ du^ 


(B.1.10) 


For a plate or a plane stress element lying in the x-y plane, 

=0. In addition, Kirchoff’s hypothesis (normals remain 

z y 21 z X 

3 V 3w 

normal to the middle surface) is equivalent to setting 3^ " “ 3^ 

Bu Bw 

and gj = " so that (B.1.9) becomes 


,, 1 ;-r.,0r 2 ,3w,2, .,0, 2 ,3w^2, 3w 3w-,,. 

= V /iN [o) +(-5—) ] + N [o) +(-5—) ] + 2N }dA, 

G 2 ^ x"- z ^3x-^ ■' y"- z '■3y^ ■' xy 3x 3y 


(B.1.11) 


where and are the membrane stress resultants acting in the 

X* y xy ® 

prebuckling state, and A is the middle surface area of the element. 


B.1.5 


In DESAP 2 the stress resultants are taken as constant within each 
element (the resultants acting at the center of the element are used); 
consequently they can be taken outside the integral sign in (B.1.11). 

The geometric stiffness matrix [G^] of an element is defined as 

Ug = ^u.}’^[G.]{u.} , (B.1.12) 


where {u^} is the nodal displacement vector of the element associated^ 
with buckling deformation. Comparing (B.1.12) with (B.1.11), we con- 
clude that the geometric stiffness matrix has the form 


1 r r (2)^x1^ r ( 3 )ixt 0 

[G.] = [g> + [g> + [g> 


(B.1.13) 


fk) 

The unit geometric stiffness matrices [g^ of the element can be 
calculated from 


. (||)2]dA , 

A 


{u.}^[gP^J£u.}= /[0)2 . (|^)2^dA , 

A 


(B.1.14) 


||dA . 


It is conventional to neglect for plate elements, but all the terms 
in (B.1.14) are included for the plane stress element. 

In the case of one-dimensional elements parallel to the x-axis, 
we use 


= 




^8w. 2-, j 

(^) ]dx 


(B.1.15) 
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where is the axial tensile force prior to buckling, and L denotes 
the length of the element. The axial force is again taken to be con- 
stant throughout L, which enables us to write 

[G.] = , (B.1.16) 

where the unit stiffness matrix is obtained from 




(B.1.17) 


Equation (B.1,15) can be derived from (B.l.ll) only by assuming that 
the element does not undergo torsion in the course of buckling. As a 
result, DESAP 2 is not capable of handling torsional or lateral- 
torsional buckling of beams. The geometric stiffness matrix for the 

latter is very complex see Ref. [5], p, 7.2-8 -- and does not readily 

lend itself to the redesign process. 

The unit geometric stiffness of each element is calculated in 
DESAP 2 as soon as the element data is read in, and is stored on an 
auxiliary storage device. Prior to each buckling analysis, the unit 
matrices are read back into the core, and the geometric stiffness 
matrices of the elements are reconstituted according to (B.1.13) 
or CB.1.6). 
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B.2 Stress Constraints 

At the present time, the only practical means of handling stress 
constraints for large structures seems to be the concept of fully 
stressed design (FSD) . In a fully stressed structure each element 
reaches its maximum permissible stress level under at least one load 
condition, unless it is governed by the minimum size constraint. FSD 
does not generally coincide with the minimum weight design, except for 
statically determinate structures, but the differences in weight are 
small in most cases. 

A typical failure criterion of an element can be expressed in the 
general form 

f({N.}, {N*},A.) = 1 . (B.2.1) 

* 

where {N^}is the intemalforce vector of the element and {N^} contains 
the allowable forces. Equation (B.2.1) may represent a criterion for 
any kind of failure, such as yielding, fracture or local instability. 

f 

Let A^ be the element size for the current design, and A^ the size 
predicted for next (improved) design. The corresponding internal forces 
are denoted by {N^} and {N^}, respectively. For the sake of clarity, 
we limit the discussion at this time to a single load condition, and 
assume that each element is subject to a single failure criterion. 

If the improved design is to be fully stressed, it must satisfy 

{N*},a’) = 1 . (B.2.2) 

The main difficulty in using (B.2.2) is that it requires a knowledge of 

I 

{N^}, i.e. the changes in the nodal forces of each element caused by the 


L 
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redesign. This information can be acquired only by an inversion of the 
structural stiffness matrix, which is impractical for large structures, 
since the banded form of the structural stiffness matrix is lost upon 
its inversion. 

The problem is greatly simplified if we assume that the nodal forces 
do not change upon redesign. Equation (B.2.2) then becomes 

f({N.}, {N*},a’) = 1 , (B.2.3) 

f 

which can immediately be solved for one element at a time. We have 
now arrived at the stress ratio method of redesign, which is adopted 
for DESAP 1. 

The assumption of no change in the nodal forces is valid only for 
statically determinate structures under dead (size -independent) loading, 
in which case a single redesign will result in FSD. For all other pro- 
blems, equation (B.2.3) is an approximation, and it must be applied 
iteratively, updating {N^} each time, before FSD is reached. If size- 
dependent loading is present (e.g., gravity or thermal loads), the 
approximation may be poor and cause difficulties of convergence. 

When several load conditions and design criteria are used, the 
element size is calculated for each combination of loading and design 

I 

formula, and the largest value is chosen as A^. 

Once the new size of the element is known, the corresponding de- 

I t 

sign variable is determined by (B.1.2): D^ = A^/n^, or the minimum 

f * * 

size constraint: D « D (D„ is the prescribed lower limit), whichever 

m m m * 

is larger. If the design variable is common to several elements. 
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i.e., i£ equal or proportional size constraints are used, the maximum 

f 

value of D is used; that is 
m 

D = max A./n^ or D = D , (B.2.4) 

m . 1 1 m m 

i£m 

whichever is greater. The notation iem indicates that maximum is ob- 
tained from all elements i that share the design variable number m. 

The ratio 


R 


m 


d’/d 

m 


is called the stress ratio of the design variable. In order to 
monitor the progress of the design sequence, the maximum and minimum 
stress ratios of the structure. 


R = max R 
max m 

m 


and 


R . * 

min 


are computed after each redesign cycle, 
be a stress-critical design if 


min R , (B.2.6) 

mm' ^ 

The new design {D } is said to 


1-6<R <1+6, (B.2.7) 

where 6 is a small parameter prescribed by the user. Similarity the 
design is considered to be fully stressed if, in addition to being 
critical, it also satisfies 


1-6<R. <1+6. (B.2.8) 

- min — V / 

The values of R , R . and the corresponding design variable and 
max' min r & * 

load numbers are printed out after each redesign cycle. 
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B.3 Buckling Constraints 

Buckling analysis is a very expensive operation that consumes 
several times more computer time than analysis of the prebuckling state. 
Moreover, the cost of buckling-constrained design is proportional to 
the number of load conditions used, since each load condition results 
in a different geometric stiffness matrix, thus requiring separate 
analysis . 

In order to keep the computer costs within reasonable limits, 

DESAP 2 allows only one load condition to be used for buckling con- 
straints, This restriction does not apply to stress constraints, where 
the number of load conditions is limited only by the available core 
storage. 

Constraints on general buckling loads are handled in DESAP 2 
by the technique described in References [3,4], The procedure is not 
only applicable to buckling loads, but can also be used for displacement 
constraints as was done in DESAP 1 [1], 

In the analysis of buckling, the loads applied to the structure 
are taken to be proportional to a single load parameter p. Denoting 
the critical values of p by p^, p^ — — * ' * ■* constraints on the 

buckling loads are 

Pp ^ P*> r = 1,2, . . . , R, (B.3.1) 

where p* is the prescribed lower bound and R equals the number of de- 
grees of freedom in the finite element model. The load vector of the 
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structure at buckling is thus given by 

{Q*} = p*{Q} , (B.3.2) 

where {Q} is the load vector used in prebuckling analysis. If {Q} 
represents the design load with respect to stress constraints, then p* 
can be interpreted as the additional factor of safety against general 
buckling. 

The constraints (B.3.1) can.be divided into two categories: if 

the equality sign governs the optimal design, the constraint is said to 
active ; if the inequality occurs, the constraint is passive . This 
division is not known, of course, until the final design is reached, 
but for notational convenience we presume that the active constraints 
are listed first, so that (B.3.1) can be replaced by 


Pr 

= p; . 

' ’ '*act 

Pr 

> P; , 

' = ^ “ 


(B.3.3) 


where R denotes the number of active buckling constraints. DESAP 2 
ac L 

allows for R < 2, i.e. it assumes that number of active buckling 
jnodes never exceeds two. Theoretically it is possible for the optimal 
design to be governed by three or more modes (all having the same buckling 
load), but such cases have not yet been encountered. 

In addition to buckling constraints, we must also account for 
limits on element sizes: 


A. > A*, 
1—1 


i.e. 


D > D* . 
m — m 


(B.3.4) 
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A passive design variable is governed by the miniiman size constraint 

* 

in the final design, i.e., D = D , whereas an active design variable 

^ mm 

* 

is determined by the displacement constraint, in which case D* > D . 

mm 

In the mathematical treatment that follows, it is convenient to 
replace (B.3.1) and (B.3.4) by the equality constraints 


P 


r 




(B.3.5) 


where a and b are to be viewed as variables free from constraints, 
r m 

The design objective is to minimize the total weight of the struc- 
ture W = EW. . Substituting for W. from (B.1.1) and (B.1.2), we have 
i ^ ^ 

W = Ep-n-D . (B.3.6J 

X m ^ 

Minimizing (B.3.6) subject to (B.3,5) is equivalent to making the 
following function stationary: 


V = Ep-n-D - EA (p -a^) - Ep (D -b^) 
1 m 1 - r^^r r'^ m 


(B.3. 7 ) 


where A^ and are non-negative Lagrangian multipliers . The operations 
3V/3D^ = 0, 3V/3a^ = 0 and 3V/3b^ = 0 yield, respectively. 


E p-n- - EA_p - y D = 0 , 
1 1 j. jr r,m m m 


A a = 0 , 
r r ’ 


y b = 0 , 
m m 


where we used the notation ( ) ni ” 


(B.3.8) 

(B.3.9) 

(B.3.10) 



Inspection of (B.3.S), (B,3,9) and (B.3 . 10)reveals that 


> 0 


X 


T 



m 



if p = 

if Pj. ^ 


★ 



if D 


m 


if D 


m 


* 



4r 

> D 

m 


(active constraints) 
(passive constraints) , 

(passive design variables) 
(active design variables) , 


which enables us to rewrite the optimality criterion (B.3, 8) as 

* 


Z p.Tl- 
1 

lem 


^ KK m ^ 

^ r r ,m 
r act ' 


>0 if D = D 

— mm 

= 0 if D > D* 

m m 


(B.3. 11) 


The notation r act shows that the sum is to be taken over the active 
displacement constraints. Introducing the unit weight of the 
design variable 


0. ■ '’i"l ■ 

lem 

(B.3. 11) becomes for the active design variables 


(B.3. 12) 


- R act 

— ZAP =1 
p , r r,m 
'^m r= 1 


(B.3. 13) 


The redesign formula for the active design variables is obtained 
directly from the optimality criterion. Multiplying both sides of 
(B.3. 13) by (l-a)D^, where a is a constant to be determined later, and 


rearranging terms, we get 
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The fimction of the design algorithm is to solve the last equation by 
successive iterations. The redesign formula used in each iterative step 
is 


D Ract 

D* = a D + (1-a) — E X p , 
m m ^ p T r^r,m 

r=l ^ 


(B.3.14) 


where is the current design variable, and represents its improved 

value. The constant a can now be recognized as the relaxation factor , 
which is utilized to control the convergence characteristics of the 
iterative procedure. 

Each iterative design, i.e. each application of (B.3.14) is followed 
by an analysis of resulting structure, so that X^ and p^ ^ can be updated. 
The computation of these parameters will be discussed next. 

Buckling is governed by the matrix characteristic value problem 


[K]{u} = p[G]{u} , 


(B.3.15) 


where [K] and [G] are the elastic and geometric stiffness matrices of 
the structure, respectively, and {u} is the nodal displacement vector 
due to buckling. The gradients of the buckling parameters can be 
shown to be [3,6] 


f [K ] - P [G 
\ ^ ,m-‘ ^r*- ,m 


]){u 


Cr) 


r,m 


{u 


Cr) 


}^[G]{u*^^^; 


CB.3,16) 


r I* 'j 

where {u^ is the, nodal displacement vector associated with p^. 

It is assumed, as was done in the stress-constrained design, that 
the changes in the internal forces between two successive redesign 
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cycles are negligible, in which case we can set [G ] = [0]. In addi- 
tion, the buckling inodes are normalized with respect to [G] , so that 
(B.3.16) becomes 

Pr,m ' . (B.3.17) 


Noting that only the elements that share the design variable con- 
tribute to [K , (B.3.17) can be rewritten as 


P 


r,m 


iem 




(B.3.18) 


Utilizing (B.1.2) and (B.1.3), the derivatives of the element stiffness 
matrices in (B.3.18) are 


[K. ] = 


n,([kf J] + n,A " [k|^^]) . 


3D 


m 


1 1 


(B.3.19) 


Since the unit stiffness matrices of each element are computed and 
stored at the beginning of the program (they do not change upon redesign), 
the gradients of the buckling load parameters are readily^ calculated once 
the buckling modes of the current design are available. 

It must be reiterated that (B.3.18) is valid only when the internal 
forces prior to buckling are independent of the design variables. This 
condition is satisfied only if the prebuckling state is statically deter- 
minate and size-dependent loads (gravity loading and thermal stresses) 
are absent. 

If the prebuckling state is statically indeterminate, (B.3.18) 
is only an approximation, which would not necessarily lead to a true 
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optimal design. The use of the equation is however, compatible with 
the stress ratio method of redesign, where the changes of internal 
forces were also neglected. Note that the use of the exact gradient 
expressions (B.3.16) is precluded for the reason stated in Sec. B.2: 
it would require the knowledge of {N^ ^}, the gradients of the internal 
force vector of each element. 

If the prebuckling state of the structure is statically determinate 
and if the loading is not size-dependent, the user of DESAP 2 should 
specify this in the input data. Since the geometric stiffness will be 
unchanged during redesign, the program will compute [G] only once, 
namely after the first prebuckling analysis. Otherwise [G] will be 
recomputed in each redesign cycle from the unit geometric stiffness 
matrices of the elements. 

The Lagrangian multipliers are chosen such as to make the improved 
design buckling-critical , i.e., they are calculated from the condition 
^r ~ ^r^ ^ ^acf ^sing the notation 6p^ " 

6Dm ” change in the nodal displacements can be estimated 

from the linear approximation 
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Z 

m=l 


6D 

,m m 
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Passive design variables are governed by the minimum size constraint 

after redesign. Hence, 6D = D* - D , making their contribution to 
® m m m 
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m pass 


p (D*-D ) 
^r,m^ m m^ 


(5p ) 

'• ^r-'pass 


(B.3.21) 
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The improved values of active elements are given by (3.3.14), which 
results in 


M 


act 


C5p,)act = ^ ^ • CB.3.22) 

m act ’ s=l ’ 


Tlie total change in is therefore. 


<5Pr = f^Pr^pass " 


f^Pr^act 


(B.3.23) 


Setting 6p^ = P* - and substituting (B.3.21) and (B.3.22) in 

(B.3.23), we obtain the following simultaneous equations for the 

Lagrangian multipliers X,s=l,2,...,R : 

s act 


act 

(1-a) E 
s=l 


M p p 
^ ^s,m^r,m ^ 


m act 


m 
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M 

E 


= (1-a) E p D - E p (D*-D ) + p* - p , 
^ ^ ^ ^r,m m ^r,m^ m m^ ^r ^r " 

m act m pass 


1,2, ..., R 


act 


The solution of (B.3.24) requires the prior knowledge of the number 

of active buckling constraints and the active-passive identities 

of the design variables. As this information is generally not available, 

the redesign is carried out by an iterative procedure diagramed in 

Fig. B.3.1. As noted previously, a maximum of two buckling modes 

are assumed to be active at any stage of the design procedure, i.e., we 

take R . < 2. 
act — 
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The flow diagram in Fig. B.3.1 is somewhat simplified, as it 
omits the paths that are followed if the equations for and X^ are 
singular or have a small determinant. In general, if X is 

indicated by the equations, the corresponding buckling constraint is 
taken as passive, i.e. is treated in the same manner as if it were 
negative . 

After each buckling analysis, but before buckling-constrained 
redesign, the buckling ratios of the current design 

Qj. = PVP^ , (B.3.25) 

th 

are computed and printed. A buckling mode (r mode) is considered to be 
potentially active only if 


Q 


r 


> wR 


max ^ 


(B.3.26) 


where 03 < 1 is a user-supplied constant and R is the maximum stress 
ratio defined in (B.2.6). If the above inequality is violated, the 
corresponding mode is ignored in the subsequent buckling-constrained 
redesign operation. 

A design is said to be buckling- critical if 


1 




ax 


max Q < 1 + 6 , 
r 


(B.3.27) 


where 6 is a small parameter that was also used in defining a stress- 
critical design see (B.2.7). The design is considered acceptable 


L 
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if one of the following conditions are met: 

1) All design variables are passive and 0 <1+6. 

Tn3.x 

2) The design is buckling-critical and the optimality 
criterion CB.3.13) is satisfied within a prescribed latitude: 

R ^ 

, act 

1 - 56 < — Z X p <1 + 56 (B.3.28) 

p T 3T r ,m — 

^m r= 1 ^ 


for all active elements. Experience with the program has shown that 
the ^^latitude" in (B.3.28) should be considerably less stringent than 
in (B.3.27); hence the use of 56 in the last inequality. 

The quantity 


, ^act 
— Z 
^m r=l 


X p 
r^r,m 


which we call the optimality index of the design variable D^, is 
printed for each design variable after every redesign cycle, together 
with its active-passive classification. 
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B.4 Stress and Buckling Constraints 

If both stress and buckling constraints are present, each 
redesign cycle is divided into two parts. First, the structure is 
redesigned with respect to the stresses only; the buckling con- 
straints are ignored in this phase of design. The next step is the 
buckling-constrained design, where the design variables just obtained 
from the stress-constrained phase are used as the minimum size constraints. 

The buckling-constrained redesign phase will be skipped if the 
analysis of the current design shows that 




< (jl3 R 
ax — max 
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where R and 0 are the maximum stress and buckling ratios de- 
max jnax 

fined in (B.2.6j and (B.3.27), respectively. The constant u) < 1 is 
also used in conjunction with (B.3.26) in choosing the potentially active 
buckling constraints. 
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B.5 Uniform Scaling Operation. 

We recall that the redesign equations for both, stress and 
buckling constraints, were based on certain simplifying approximations. 

The stress ratio redesign formula (B.2.3) and the gradients of the 
buckling parameters (B.3.18) were both derived under the assumption that 
the internal forces are unchanged during a redesign cycle. In addition, 
the predicted changes in the buckling parameters (B.3.20) were 
linearized. As a consequence of these approximations, the redesign 
process becomes an iterative procedure consisting of repeated applications 

of the redesign formulas until convergence is achieved. 

It is frequently advantageous to interrujit the design procedure by 

the so-called uniform scaling operation,, where all the design variables 
are changed by the same scale factor y: 

D • » y . (B.S.l) 

m n 

The scale factor is calculated from the condition that the scaled design 
{D*} should be critical, i.e. 

1-6 ^ max(R^^, ^1 + 6. (B.5. 2) 

DESAP 2 gives the user the option of using the uniform scaling 
operation whenever the current design is not critical. Once the design 
has been scaled to the critical state (more than one scaling operation 
may be required to achieve this), the redesign equations will be applied 
in the usual manner. 

This method of design offers two advantages over the use of the 
redesign equations above. Firstly, the usq of the scaling operation 
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results in a sequence of critical designs which are very useful in 
monitoring the design process. In particular, the weight comparison of 
the critical designs can be used to terminate the design whenever the 
weight reduction becomes small, or ceases altogether- This is especially 
important when approximate optimization techniques, such as the fully 
stressed design principle, are used. 

The second advantage of using scaling is that it prevents inter- 
mediate designs from departing excessively from the critical state. 

This in turn has a stabilizing influence on the convergence of the 
design process, and in some problems it even makes the difference be- 
tween a convergent process and no convergence at all. 

The scale factor for the stress constraints is simply 

, (B.5.3) 

where R is the maximum stress ratio obtained by the stress ratio 

method -- see (B.2.5), The scale factor is exact that is, the 

resulting design will be precisely stress-critical if the internal 

forces remain unchanged upon scaling. This condition is satisfied, 
apart from statically determinate structures, if all element stiffness 
matrices have the form 

[K.] = [k.]A^ , (B.5.4) 

where n is common to all elements, and if the loading is size-independent. 

It can easily be shown that if (B.5.4) is satisfied, the 
buckling -constrained uniform scaling is also exact, the scale factor 
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being 




(B.5.5) 


where is the maximum buckling ratio defined in CB.3.27). 

For the problems where the scaling operation is inexact , the 
scale factor for the buckling constraints must be computed from (B.3.20). To 
obtain the scale factor ^ which would make the r buckling load 
critical, we substitute 6p^ = p* - p and (^D = D' - D = (y^^'^-l)D , 

^ ^ *■ ^ |fl TTI ^ ^ ^ TTl ^ 


m 


m 


m 


obtaining 


.(B) 
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P* - P = (y^ -1) 2 p D 

r r , ^r,m ra 

m=i 


Solving for the scale factor, we get 
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yt®) = (p*-p )/( Z p D ) + 1 

•^r j. f j.-'/ /r,m m' 

m=i 


(B.5.6) 


A scale factor is computed for each potentially active buckling 
constraint from (B.5.6), and the maximum value is used for y^^^, 


i.e. 


(B) (B) 

y'- •' = max y^ , 

r 


r = 1, . . .,R 


act 


(B.5.7) 



The factor used in the scaling operation (B.5.1) is the larger of 

j (B) • 

and y ' , i . e . 


y = max(y^^\y^®b 


(B.5.8) 


If the structure meets the conditions for exact scaling, the user 
should specify this in the input to DESAP 2 , together with the inertia 
exponent n in (B.5.4). The buckling scale factor will then be 
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computed from (B.5.5) and the iuialysis of the scaled structure will 
be omitted, resulting in a considerable saving in computer time. 
Otherwise, CB.5.6) and (B.5.7) are used to calculate and the 

scaled design will be reanalyzed and if not critical, will be scaled 
again. 

DESAP 2 gives the user the option of dispensing with scaling 
altogether. The no-scaling option could be used if the scaling oper- 
ation is inexact and the intermediate critical designs are of no in- 
terest. This could save a substantial amount of computer time, but at 
the expense of losing some control over the convergence of the design 
procedure. 

Scaling should not be used if the loading is dominated by thermal 
or gravity loads and the element stiffness matrices have the form 
[K^] = [k^]A^ (linear size-stiffness relationship). It is readily seen 
that under these circumstances a uniform scaling operation is entirely 
ineffective, since it leaves the stresses and buckling loads unchanged. 

If the uniform scaling option is chosen, the weight of each critical 

design is calculated and the smallest of these weights, W . , is stored. 

min 

The design procedure is terminated whenever 


(W-W . )/W . > e 

^ min*^ min — 


(B.5.9) 


where W is the weight of the current critical design and e is a small 
prescribed constant. This cut-off criterion prevents further redesign 
operations if these are going to result in a weight increase. 
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The scaling operation is always by-passed if y > 2, in order to 
avoid the corresponding large weight increase. Experience has shown 
that by preceding uniform scaling with a redesign operation, a lower 
weight is obtained than from the scaling-redesign sequence. It should 
be pointed out that a large scale factor would arise only during the 
first design cycle if the initial design is poorly chosen. 
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C. ORGANIZATION OF THE PROGRAM 

C. 1 Overlay Tree 

The structure of the program is illustrated with the overlay tree 
in Fig. C.1.1 Details of the element subroutines are shown separately 
in Figs. C.1.2 to C.1.7. 

Additional elements can be added to the program by replacing the 
call NOELEM in ELTYPE with the new element subroutines. The element code 
numbers (MTYPE) 5 and 8 have been reserved for this purpose. 




Figure C. 1 . 1 

Overlay Tree for DESAP 2. 
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Figure C.1.2 

Overlay Tree for Bar Element 


Figure C.1.3 

Overlay Tree for Beam Element 
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Figure C.1.4 

Overlay Tree for Plane Stress Element 
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Hgure C.1.5 

Overlay Tree for Shear Panel 
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Figure C . 1 .6 

Overlay Tree for PI ate -She 11 Element. 
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Figure C.1.7 

Overlay Tree for Boundary Element 
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C.2 Element Subroutines 

The element subroutines of DESAP 2 are constructed in such a 
manner as to facilitate modifications by the user. Each element sub- 
routine package consists of two parts: element data processing, and 

stress recovery combined with stress ratio redesign. The element data, 
including the design information, is read in by the processor sub- 
routines shown by the dashed lines in Figs. C.1.2 to C.1.6. Provision 
has been made for more than one such processor subroutine with each 
element package, the different processors being distinguished by their 
construction code numbers (KODE) , This scheme allows the user to pro- 
vide his own processor subroutines which will best satisfy his special 
design requirements. 

After the data for an element is read in, it is immediately pro- 
cessed and the results placed on auxiliary storage devices. Standard 
computations, such as the formation of the unit elastic and geometric 
stiffness matrices, transformation of coordinates, etc , are handled 
by calling the appropriate element computational subroutines , identified 
by the solid lines. No modification by the user is required here, since 
these sub-routines are independent of the construction details and de- 
sign criteria. 

The processed element data consists of: 

1) Element unit elastic stiffness matrices and load vectors. 

2) Element unit geometric stiffness matrices. 

3) Element unit force recovery matrices and force vectors. 

4) Data associated with failure and local instability criteria. 

The first three sets of data are stored in a standard form (common to all 
elements), whereas the format of the failure and local instability data 
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I 
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is completely flexible it may, for example, contain parameters of the 

failure equations (as in the subroutines presently provided), or purely 
numerical data in tabular form. 

Following each prebuckling analysis of the structure, the element 
nodal forces are recovered by subroutine STRSC. The appropriate design 
subroutine is then used to compute the stresses and carry out stress 
ratio redesign. The design subroutines, also identified by the dashed 
lines in Figs. C.1.2 to C.1.6, must be compatible with the corresponding 

processor subroutines that is, if the user supplies his own processor, 

he must also have a matching design subroutine, both labelled with the 
same construction code. 

Note that the boundary element. Fig. C.1.7, is not subject to 
redesign and, therefore, lacks a construction code and design sub- 
routine. 

The user-supplied processor and design subroutines are to be in- 
serted in the place of NOELEM and RETURN, respectively, shown in 
Figs, C.1.2 to C.1.6. 
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C. 3 Storage Requirements 

The common in-core storage blocks used in DESAP 2 are: 

1) COMMON/JUNK/ 368 locations is used for storage of miscellaneous 

data, as its name implies. 

2) COMMON/ELPAR/ 33 locations serves for storage of parameters that 

control the execution of the program. 

3) COMMON/UNITS/ 12 locations contains the assignments (numbers) of 

the input-output units. 

4) COMMON/EM/ 5548 locations is used primarily for processing of 

element data. 

5) COMMON/CONTR/ 29 locations contains control data for the redesign 

operations . 

6) Unlabelled storage area A(n). This is the main working area of the 
program and its dimension largely controls the allowable size of the 
structure. The capacity of the program can be adjusted by changing 
^^n" in the following statements at the beginning of the main program: 

DIMENSION A(n) 

REAL* 8 AD (n/2) 

EQUIVALENCE (A(l), AD(1)) 

MTOT = n 

The availability of the required storage in A(n) is checked at various 
stages of the program, and if insufficient, an error message is printed 
(subroutine ERROR) and the execution of the program terminated. Tlie 
minimum value of ^’n" required in each subroutine where A(n) is used, is 
listed in Table C.3.3. It is advisable, however, to use the largest 
possible value for ”n", since this would reduce the time spent on the 
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: :\-msfer of data between the auxiliary storage devices and the core. 

k.:OMMON/COMPL 852 locations is used only in the plate-shell element 

subroutines. 

The usage of the common blocks in the various subroutines is shown 
in Table C.3.4. 


L 
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Parameter 

Description 

LL 

Number of load conditions subject to stress constraints 

M 

Max (,LL,2') 

MM 

2*NEQB+ CMBAND-1)/NEQB 

MBAND 

Band width of structural stiffness matrix 

NELTYP 

Number of element types in the structure 

NEQ 

Number of degrees of freedom for the structure, i.e., 
number of equations 

NEQB 

Number of equations in a block 

= min{(n-4*LL)/[4*(MBAND+LL)],n/[4*(MBAND+M) + l] ,NEq} 

ND 

Number of degrees of freedom for an element (see Table C.3.2) 

NG 

Number of unit geometric stiffness matrices for an element 
(see Table C.3.2) 

NI 

Dimension of element design information array 
(see Table C.3.2) 

NS 

Number of stresses calculated for an element (see Table C.3.2) 

NU 

Number of unit stiffness matrices for an element 
(see Table C.3.2) 

NUMDV 

Number of design variables 

NUMEL 

Number of elements in structure 

NUMNP 

Number of nodal points in structure 

NUMGEO 

Number of different geometric properties 
used for an element type 

NUMTC 

Maximum number of temperatures for which material 
properties of an element type are specified 

NUMFX 

Number of fixed end-force sets (beam only) 

NV 

Number of unit load vectors for an element 
(see Table C.3.2) 

NW 

Number of initial stress vectors for an element 
(see Table C.3.2) 


Table C.3.1 


List of Parameters that Determine the Storage Requirements . 







Parameters 


Element 

NU 

NV 

NW 

ND 

NS 

NI 

Truss 

1 

1 

1 

6 

1 

5 

Beam 

2 

2 

2 

24 

12 

11 

Plane 

1 

2 

1 

12 

15 

10 

Shear 

1 

1 

1 

12 

4 

3 

Plate-Shell 

2 

2 

1 

24 

6 

5 

Boundary 

1 

1 

1. 

6 

2 

0 


Table C.3.2 


Values of Storage Parameters for Various Elements 
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Subroutine 

Minimum Required Value of "n" 

INL 

6*NUMNP+6*LL+2*NEQB*LL 

ADDSTF 

4* (MBAND+LL)*NEQB+4*LL 

USOL 

4* (MBAND+M) *NEQB+NEQB 

PRINTD 

6*NUMNP+6*M+2*NEQB*M 

STRESS 

3*NUMDV+4*LL+2*NEQB+NEQ 

MULBAN 

NEQB* (MBAND+LL) + 2*MM*iWEC 

DERV 

2*NUMDV+NEQB 

BDESIN 

6*NUMDV+4*LL 

ELTYPE 

NUMDV+10*NUMNP+m, where depends on the element type 

[see below) 


Subroutine 

Minimum Required Value of ”m^' 

TRUSS 

C2+5*NUMTC) *NUMMAT+2*NUMGE0 

BEAM 

10*NUMGEO+6*NUMMAT+12*NUMFX 

PLANE 

(2+8*NUMTC) *NUMMAT+5*NUMGE0 

SHEAR 

C2+4*NUMTC)*NUMMAT 

SHELL 

(2+7*NUMTC)*NUMMAT 


Table C.3.3 


Jll 


Core Storage Requirements for A(n) 

The storage parameters are defined in Table C.3.1, 
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MAIN 

INL 

CALBAN 

ELSTIF 

ADDSTF 

STRESS 

STRSC 

BANAL 

ADGSTF 

DERV 

DESIGN 

BDESIN 


ERROR 

TRUSS 

RUSS 

FTRUSS 


DTRUSS 

TGEOM 

BEAM 

TEAM 


NEWBM 

SLAVE 

DBEAM 

BGEOM 


PLANE 

ELAW 

QUAD 

FORMB 


PLNAXl 

DPLANl 

PLNAX2 

DPLAN2 


SHEAR 

PANEL 

FPANEL 

DPANEL 



SLST 

SLCCT 

DSHELl 

SHELGl 


SHELG2 

BOUND 

CLAMP 



Table C.3.4 


Usage of Common Blocks in Various Subroutines. 
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C.4 Requirements for Auxiliary Storage Devices 

DESAP 2 requires nine sequentially accessible, auxiliary storage 
units* The units, together with the required storage capacities, are 
listed below. The storage parameters can be found in Table C.3,1. 

11 is used for permanent storage of element load multipliers and 

the minimum allowable design variables. 

No. of locations = 4* (LL+NUMDV) 

12 is used for: 

a) temporary storage of element stiffness matrices, 

b) scratch file during solution of simultaneous equations, 

c) storage of displacement vectors, 

d) storage of the coordinate vectors during buckling analysis, 

e) temporary storage of the derivatives of the buckling parameters. 
No. of locations = the larger of the following: 

a) Z [2+ND*(ND*2+9)] , 

NUMEL 

b) 2 * |mband/neqb|*mband*neqb, 

c) 2*NEQ*M, 

d) 2*NUMDV, 

13 is used for; 

a) scratch file during solution of simultaneous equations, 

b) temporary storage of element geometric stiffness matrices, 

c) storage of structural stiffness matrix and the coordinate 
vectors in the solution of equations during buckling analysis, 

d) storage of mode shapes during buckling analysis. 
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II 


No. of locations = the larger of the following: 

a) 2*NEQ* (MBAND+LL) , 

b) E [2+ND*(l+2*ND) , 

NUMEL 

c) 2*NEQ*(MBAND+2), 

d) 4*NEQ. 

18 is used for permanent storage of the equation number matrix, 

control parameters for element subroutines (NPAR) , element unit 
force recovery matrices, element stress -constrained design in- 
formation, and unit weights of the design variables. 

No. of locations = 6*NUMNP+14*NELTYP+NIJMDV 

+ E [7+ND+NI+(ND*NU+4*NW)*2*NS] . 

NUMEL 


19 is used for: 

a) scratch file in the assembly of structural stiffness and 
geometric stiffness matrices, 

b) scratch file in the solution of simultaneous equations (during 
prebuckling and buckling analysis), 

c) scratch file during buckling analysis. 

No. of locations = the larger of the following: 

a) 1 (NUMEL)^"^^ |*[ND* CND*2+9)+2] 

b) |MBAND/NEQB|*2*MBAND*NEQB, 

c) 4*NEQ. 
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no is used for storage of structural stiffness matrix, load vectors 

and geometric stiffness matrix (statically indeterminate structures 
only) . 

No. of locations = 2*NEQ* (2*MBAND+LL+2) . 

111 is used for: 

a) permanent (indeterminate structures) or temporary (determinate 
structures) storage of element unit geometric stiffness matrices, 

b) permanent storage of the structural geometric stiffness matrix 
(statically determinate structures only). 

No. of locations = the larger of the following: 

a) Z (4+ND*ND*NG) , 

NUMEL 

b) 2*NEQ*MBAND. 

112 is used for permanent storage of element unit stiffness matrices 

and load vectors, and permanent storage of structural load vectors 

due to size-independent loads. 

No. of locations = Z [7+ND+2*ND* (ND+NU+4*NV) ] 

NUMEL 

+ 2*LL*NEQ. 

113 is used for: 

a) storage of coordinate vectors used in buckling analysis, 

b) scratch file during buckling analysis. 

No. of locations = the larger of the following: 

a) 4*NEQ, 

b) 2*NEQ*(MBAND+2) . 
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The assignments of the auxiliary storage devices, together with 
those of the input-output units (IR = card reader, IW = printer, and 
IP = card punch) , can be changed by altering the assignment numbers at 
the beginning of the Main Program. 
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We denote the load vector of the structure representing the 
5/th load condition by {Q^ and the corresponding nodal displacement 
vector by {u^^^}. The equilibrium equations governing the prebuckling 
state thus are 

[K][U] = [Q] , (C.5.1) 

where [K] is the stiffness matrix of the structure, 

ru] = ... , 

(C.5.2) 

[Q] = ... , 

and L represents the number of load conditions. 

The assembly of the equations and their solution, carried out by 
ADDSTF and USOL, respectively, are essentially the same as in the SOLID 
SAP program. The only significant difference is that DESAP 2 uses 
double precision arithmetic, whereas single precision was employed in 
SOLID SAP. Reconversion into single precision is not recommended, even 
if your machine carries above-average number of digits, because roundoff 
errors in the analysis tend to be magnified during the redesign cycle. 
The cumulative buildup of error may even lead to a non-convergence of 
the design process in large structures. 

The equations, i.e., the matrices [K] and [Q] , are formed in 
blocks in the unlabelled common area A(n), and stored on auxiliary 
storage devices. The banded structure of the stiffness matrix is 
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[ 


exploited throughout the formation of the equations and their solution. 
The size of a block is determined automatically by the computer once 
the bandwidth and the number of load vectors have been established* 
There must be sufficient space in A(n) for at least two equations. 

The equations are solved by the Gaussian elimination procedure, 
one block at a time. The structural stiffness matrix is not destroyed 
during the solution procedure. 
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C.6 Solution of Buckling Equations 

As pointed out in (B.5.15), buckling analysis leads to a matrix 
eigenvalue problem of order R (R = number of degrees of freedom) : 

[K]{u} = p[G]{u} . (C.6.1) 


DESAP 2 uses an iterative Rayleight-Ritz method [7] to extract only the 


potentially active eigenvalues (load parameters) p^^^, 


»act’ 


and the corresponding eigenvectors (buckling modes) {u^^^}, 
{u }. 


The first step is to introduce the R-dimensional coordinate vectors 
{z^ r = 1,2, . . . and to express the nodal displacement vector 
in the form 


{u} 




} . 


(C.6.2) 


where y^ are called the Rayleigh-Ritz coordinates. Equation (C.6.2) 
can also be written as 

{u} = [z]{y> , (C.6. 3) 

where [Z] = ... 

The coordinate vectors can be interpreted as the first approximation 
for the active buckling modes, whereas {u} represents an improved 
approximation. 

Once the coordinate vectors have been chosen, the Rayleigh-Ritz 
coordinates can be obtained by minimizing the potential energy of 


C.6.2 


the buckling deformations: 


E = i({u}'^[K]{u} - p{u}'^[G]{u}) . 


(C.6.4) 


Substituting (C.6.3) in (C.6.4) and applying the operations 

9E/9y =0, results in the reduced (order R ) eigenvalue problem 

s act 


[K']{y} = p[G']{y} . 


(C.6.5) 


where 

[K>] = [Z]'^[K][Z], [C] = [Z]'^[G][Z] . (C.6.6) 

Since R < 2 in DESAP 2, the solution of (C.6.5) can easily be written 
act 

down once [K*] and [G*] have been formed. 

The improved estimate for the buckling modes can now be obtained 
by substituting the eigenvectors of (C.6.5) {y^^^}, 

(R .) 

{y } in (C.6.3). Thus 

= [Z]{y‘^''^}, r = , 


or 


[U] = [Z][Y] . 


CC.6.7) 


where 

[U] . {u'^b ... 

[Y] = ... . 


and 
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It is shown in [7] that the knowledge of [U] enables us to obtain 
a better approximation to the coordinate vectors : 

[Z] = [K"^][G][U] , (C.6.8) 

which can be used in turn to obtain closer estimates of the critical 
load parameters and the buckling modes. The process is repeated until 
the change in the eigenvalues becomes sufficiently small: 

max |Ap^^Vp^^^| £ 0.005 , (C.6.9) 

r 

fr') fr'l 

where Ap^ ^ is the change in p between the last two iterations. 

Since the inversion of a banded matrix is a very inefficient 
operation, the improved coordinate vectors are not calculated directly 
from (C.6.8), but are obtained from the solution of the simultaneous 
equations 


[K][Z] = [C] , 


(C.6.10) 


where [C] = [G] [U] . 

The initial choice of the coordinate vectors [Z] can be made in 
two ways : 

(1) the coordinates are chosen by the user and punched on data 
cards in the order of the equation numbers, or 


(2) the coordinates are assigned by the computer by means of a 
random number generator in subroutine INPUTZ. 
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These coordinate vectors are used only to start the iterative analysis 
of the initial design. In the analysis of all subsequent designs, the 
initial coordinate vectors are taken as the buckling modes of the pre- 
vious design. Since the buckling modes of two successive designs do not 
change much, the number of iterations required for each buckling analysis 
will be small, with the possible exception of the initial design. 

The geometric stiffness matrix of the structure [G] is assembled 
by the subroutine ADGSTF. It differs from ADDSTF, which is used for 
the assembly of the elastic stiffness matrix [K] , in minor details only. 
There is no need to recompute [K] for the buckling analysis, since it is 
already available from auxiliary storage where it was placed during 
analysis of the prebuckling state. The simultaneous equations (C.6.10) 
are solved by the same subroutine that was employed for the equilibrium 
equations of the prebuckling state, namely USOL. A flow diagram of 
the buckling analysis algorithm appears in Fig. C.6.1. 

















C.7,1 


C. 7 Cutoff Criteria and Restart Option 

The design process is terminated whenever one of the following cutoff 
criteria are satisfied: 

1) The number of design cycles equals a prescribed number NCYCL. 

A redesign cycle is defined as an application of the redesign equations; 
uniform scaling operation, if used, does not count as a redesign cycle. 

This feature enables the user' to evaluate the progress of the design 
procedure after a few design cycles and take corrective action if neces- 
sary (particularly, adjust the relaxation parameter a). 

2) The number of successive uniform scaling operations equals a 
prescribed number NSCALE. Program termination here indicates that uni- 
form scaling is an ineffective operation and should not be used. 

3) The structural weight of critical design begins to increase, i.e.. 


(W - W . )/W . > e , 

^ min^ min — ' 


(C.7.1) 


where the terms are defined in (B.5.9). 

4) The design is acceptable, i.e., convergence is complete. An 
acceptable design satisfies one of the following optimality conditions: 
(i) The design is fully stressed and the buckling con- 
straints are not violated, i.e.. 


1-5<R <1+6, 1 

— max — 


6 


R . < 

min 


+ 6 


(C.7.2) 


Qmax 1 ^ ^ * 

where R , R • and 0 are the stress and buckling 
max* min ^ax 

ratios defined in (B.2.6) and (B.3.27), respectively, and 
6 is a prescribed (small) constant. 
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(ii) The design is buck ling -critical, the optimality 

criterion is satisfied for the buckling constraints, 
and stress constraints are not violated, i.e.. 


1 - 6 < < 1 + 6 ' 
j act 

1 - 56 < — Z X P < 1+56 for all active D , 

-Pm r=l 

(C.7.3) 

and R <1 + 6 
max — 


At the user's option, DESAP 2 will produce a restart deck just prior 
to the termination of the program. The restart deck contains the last 
values of the design variables together with their minimum allowable 
values. It can be used to replace the original design variable data 
input deck if the user decides to resubmit the program and use the last 
design as the starting point. 
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D. DESCRIPTION OF INPUT DATA 


D. 1 General Input Data 
Heading Card (20A4) 


1 


80 


HED 



HED = Any alphameric statement that is to appear as the first line 
of output. 

Program Control Card (415) 


1 

6 

11 

16 20 

NUMNP 

NELTYP 

LL 

NUMDV 


NUMNP = Number of nodal points in the structure. Include special 

nodes used for beam elements (slave nodes and points used to 
define principal axes) and boundary elements (points used to 
define directions of the elements). 

NELTYP = Number of different element types used in the structure. Count 
each construction code of an element as a separate element type. 

LL = Number of separate loading conditions for which the structure 
is to be designed. Note that only one of these loading con- 
ditions may be subjected to buckling constraints. 

NUMDV = Number of independent design variables in the structure. 


III. Design Control Card (315, 2F10.0, 315) 


1 

6 

11 

16 

26 

36 

41 

46 50 

NCYCL 

NSCALE 

KSCALE 

DELTA 

EPSIL 

KPUNCH 

KPRINT 

LBUCK 


NCYCL = Maximum allowable number of redesign cycles (see Sec. C.7). If 
NCYCL = 0, analysis of initial design only will be carried out. 

NSCALE = Maximum allowable number of successive scaling operations (see 
Sec. C.7). If blank, computer sets NSCALE = 3. 
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KSCALE = Code for uniform scaling operation. 

KSCALE < 0: omit uniform scaling operation. 

KSCALE = 0: carry^ out uniform scaling whenever the design is 

not critical, and analyse the scaled structure. Used when 
scaling is not exact (see Sec. B.5). 

KSCALE = n > 0, where n is the exponent of A^ in [Ki] = [ki]A^. 
Uniform scaling will be carried out whenever the design is not 
critical, but analysis of the scaled structure will be skipped. 

To be used only when scaling is exact (see Sec, B.5). 

DELTA = The small parameter 6 that specifies the latitude of the cut-off 
criteria listed in Sec. C.7. If blank, computer sets DELTA=0.05. 

EPSIL = The small parameter e that specifies the allowable weight in- 
crease of the structure in the weight-cutoff criterion in Sec. C.7. 
If blank, computer sets EPSIL = 0.1. 

KPUNCH = Code for restart card deck. 

KPUNCH ^ 0: punch design variable data cards for the last 

design before program termination occurs (see Sec. C.7). 

KPUNCH = 0: omit punching. 

KPRINT = Code for printout of nodal displacements, incl . buckling j^odes. 

KPRINT ^ 0: print nodal displacements after each analysis of 

the structure, 

KPRINT = 0: omit printout of nodal displacements. 

LBUCK = Code for constraints on buckling loads. 

LBUCK =0: no buckling constraints exist. 

LBUCK = n > 0: the nth load condition is subjected to constraints 

on the buckling loads. 

IV. Nodal Point Data (715, 3F10.0, 15, FIO.O) 

One card for each node; cards do not have to be in node number sequence, 
but the last card in the deck must be for the last node. 


1 

6 

11 

16 

21 

26 

31 

36 

46 

56 

66 

71 80 

N 

u 

1 U 

1 u 

1 S 

1 e 

i 9 

X 

1 Y 

1 Z 

KN 

T 


x 

' y 

' z 

' X 

' y 



1 
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N = Node number. 

ID = Motion code for nodal displacements (u^, Uy, U 2 ) and rotations 
9^, 6^, 0^) with respect to the global coordinate axes. 

ID(i) = 0: motion is permitted. 

ID(i) = 1: motion is not permitted (the corresponding degrees 

of freedom are eliminated), 

ID(i) = m, m > 1: node N is connected to a master node m. 

(used for beam elements only see Ch. E) . 

X,y,Z = Global coordinates of the node. 
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KN = Node number increment used in automatic generation of nodal 
points (see below) . 

KN = 0: automatic generation is not used. 

KN > 0: use automatic generation. 

T = Nodal point temperature. 

Note: temperature distribution in the structure is prescribed 

by nodal temperatures. The reference temperature (i.e. temper- 
ature of the stress-free state) is specified on the element cards. 

Automatic node generation if a series of nodes exists that satisfies 

the following requirements: 

(a) the nodes are equally spaced along a straight line; 

(b) the temperature distribution along the line of nodes is linear 

(c) the motion codes for all the nodes are identical; 

(d) the node numbers form the sequence n, n+KN, n+2*KN,..., 
then only cards for the first and last node of the series are required. 

The data for the intermediate nodes will be generated automatically. The 
motion code for the series will be taken from the first card ^ the node 
increment KN from the last card of the series. 

Automatic motion code generation if a particular degree of freedom is 

to be eliminated from a series of nodes, this can be achieved by using 
-1 as the motion code on the first card, and 1 on the last card of the 
series. The computer will automatically use 1 as the motion code for 
the first and all intermediate node cards of the series. 

V . Element Data 

This data differs for various element types and construction codes; it is 
described from Ch. E onwards. 

VI. Structural Load Multipliers (4F10.0) 

One card for each separate load condition; the cards must be in ascending 
order of load condition numbers. 


1 n 21 31 40 

A I B I C I D 
- <. - STR(4) — > - 


STR(l) = Fraction of element load A which is to be added to the load 
condition (structural load multiplier for element load A). 
STR(2) = Structural load multiplier for element load B. 

STR(3) = Structural load multiplier for element load C. 

STR(4) = Structural load multiplier for element load D. 
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Element loads A, B, C and D are defined with description of element 
data (Ch. E onwards). The element loads represent thermal and gravity 
loading or distributed loads of prescribed magnitude. The structural 
load multipliers enable the user to choose any fraction of the element 
loads for any load condition. 

VII. Buckling Control Card (FIO.O, 415, 2F10.0) 

This card is to be omitted if no buckling constraints are imposed 
(LBUCK = 0). 


1 

11 

16 

21 

26 

31 

41 50 

COEFFT 

MODEIN 

NMODE 

INDET 

NVEC 

ALFA 

OMEGA 


COEFFT = The constraint p* on the critical load parameters, defined in 
eqs. (B.3.1) and (B.3.2). It represents the desired factor of 
safety against general buckling for the prescribed loading. If 
blank, computer sets COEFFT = 1.0. 

MODEIN = Code for generation of the initial coordinate vectors (mode 
shapes) see Sec. C.6.1. 

MODEIN = 0: coordinate vectors are to be generated by the 

computer as random numbers . 

MODEIN = 1: coordinate vectors are read in with input data. 

NMODE = Number of buckling modes on which the buckling constraint is to 
be imposed. 

NMODE = 1: the optimal design is expected to be determined 

by a single (fundamental) buckling mode. 

NMODE = 2: the optimal design is expected to be determined by 

two buckling modes simultaneously (i.e., the lowest two buckling 
loads are expected to coalesce). If blank, computer sets 
NMODE = 1; if NMODE >- 2, computer sets NMODE = 2. 

Note: The use of NMODE = 1 is recommended for the first run, 
since it works for a vast majority of problems and requires 
somewhat less computer time than NMODE = 2. 

INDET = Code for static determinacy of the prebuckling state. 

INDET = 0: the internal forces are statically indeterminate 

prior to buckling. 

INDET = 1: the internal forces are statically determinate 

prior to buckling. 

Note: the code applies only to the load condition that is 

subjected to the buckling constraints. If INDET = 1, the 
geometric stiffness matrix of the structure is assembled only 
once; if INDET = 0, the matrix is reassembled prior to each 
buckling analysis see Sec. B.3.9. 
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NVEC = Number of coordinate vectors {z'* used in the iterative 
buckling analysis. If NVEC < NMODE or NVEC > 2, computer 
sets NVEC = 2. 

Note: the number of coordinate vectors used must equal at 

least the number of eigenvalues (NMODE) to be extracted. 

It is advisable, however, to use NVEC = 2 (the maximum allowable 
value) in all cases (even if NMODE = 1), since this results in 
faster convergence. 

ALFA = The relaxation factor a < 1 used in the buckling-constrained 
redesign formula (B.3.13). 

Note: the chances that the iterative design procedure con- 

verges uniformly increase with increasing values of ALFA. 
Consequently, it is better to use an ALFA that is too large 
(resulting in under-relaxation), than one that is too small 
(over-relaxation) . Large values of ALFA, however, may require 
a larger number of design iterations before the cutoff criteria 
are activated. In order to minimize computer time, it is ad- 
visable to start with a '^normal" value, which can be estimated 
from [3] a = n/(n+l), where n is the inertia exponent that 
dominates the structural stiffness matrix. After a few re- 
design cycles ALFA may be increased or decreased, depending 
on the convergence characteristics of the problem. 


OMEGA = The comparison parameter o) < 1 used in selecting the potentially 
active buckling constraints (see Sec. B.4). If blank, computer 
sets OMEGA = 0.8, 

VIII. Nodal Load Data (215, 6F10.4) 

One card per load condition for each node that carries concentrated loads 
or moments. Cards must be arranged in an ascending order of node 
numbers. The data must end with a blank card . If no nodal loads are 
present, use the blank card only. 


1 6 11 21 31 41 51 61 70 


N 

L 

R 1 R 

x ' y 

R 1 M 1 M 1 M 

z ' X ' y z 


R(6) 


N = Node number. 

L = Load condition number. 
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R = Concentrated forces (R ,R ,R ) and moments (M ,M ,M ) in the 

X y z X y z 

global coordinate directions. Positive values denote loads 

in the positive coordinate directions, negative values in the 

negative coordinate directions (see Fig. D.1.1). 



Note: the right-hand screw 

rule is used to determine 
positive directions for the 
moments, as shown in Fig, D.1.1. 
The same sign convention 
is employed for nodal dis- 
placements and rotations in the 
printout of the analysis. 


Fig. D.1.1 

Positive Directions of 
Nodal Forces and Moments. 


IX. Initial Coordinate Vectors (8F10.5) 

C I* 'I 

One set of cards for each coordinate vector {z^ used in buckling 
analysis. The number of vectors must equal NVEC specified on the 
Buckling Control Card. The components of each coordinate vector are 
to be read in by blocks in the order of their equation numbers? If 
the initial coordinate vectors are to be generated by the computer 
(specified by MODEIN = 0 on the Buckling Control Card), omit all cards. 

Design Variable Data 

One card for each design variable; cards must be in ascending order of 
design variable numbers. 


1 6 16 25 


N 


AOLD 


AMIN 


N = Design variable number. 

*Note that prior knowledge of the equation numbers and the number of 
equations in each block are required. This information is available 
only after the first run of the problem has been completed. 
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AOLD = Initial value of design variable. If AOLD < AMIN, computer 
sets AOLD = AMIN. 

AMIN = Minimum allowable value of design variable. 

Automatic data generation if there exists a sequence of design variables 

N = n, n+1, n+2, ... that have identical values of AOLD and AMIN, it is 
sufficient to include only the last card of sequence in the deck. The 
omitted data will be generated automatically. 
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E. BAR ELEMENT 


A typical bar element is 
shown in Fig. E.1.1. The 
element coordinate axes y and 
z coincide with the principal 
axes of the cross section, 
whereas x is the centroidal 
axis . 

The displacements are 
assumed to be linear in x; 
consequently, the axial force 
P will be constant throughout 
the length of the element. 

Since bar elements do not carry bending or torsion, the nodal rota- 
tions must be suppressed on the Nodal Point Data cards, utilizing the 
motion code ID(6) see Sec. D.l. For the same reason, only zero con- 

centrated moments are permitted on the Nodal Load Data cards. 

The basic element loads consist of gravity loading (in the three 
global coordinate directions), and temperature increase. The temperature 
increase of an element is assumed to be uniform, and is calculated from 

AT = 1/2(T.+T.) - T^^^ , (E.1.1) 

where T^ and T^ are the nodal temperatures (see Nodal Point Data), and 
T^^^ is the reference temperature specified on the element cards. The ele- 
ment load vector due to gravity is obtained by assigning half the 
gravitational force acting on the element to each node. 


E. 1 General Information 

x 



Fig. E.1.1 

TN'pical Bar Element. 



E.1.2 


The program can handle temperature -dependent material properties, 
in which case the properties should be listed for two or more tempera- 
tures, The lowest and highest of these temperatures must cover the range of 
temperatures listed on Nodal Point Data cards. The material properties 
of each element are then obtained from the listing by linear interpolation. 

The output from the analysis will consist of the axial force P for 
each element and each load condition. 

The element data deck for each construction code used must start 

with: 

IV A. Element Control Card (615) 


1 

6 

11 

16 

21 

26 30 

MTYPE 

NUME 

NUMMAT 

NUMGEO 

KODE 

NUMTC 

\TriAn/'z:'\ — 1 







NPAR(l) = Code for element type (MTYPE) . For bar elements use NPAR(l) = 1. 

NPAR(2) = Number of bar elements with the specified construction code (NUME) . 

NPAR(3) = Number of different materials used for the specified construction 
code (NUMMAT) . 

NPAR(4) = Number of different geometric (cross-sectional) properties used 
for the specified construction code (NUMGEO) . 

NPAR(5) = Construction code (KODE) . 

NPAR(6) = Maximum number of temperatures for which material properties are 
given (NUMTC) . If blank, computer sets NPAR(6) = 1. 

The remaining element data depends on the construction code used, and is 

described separately below. 
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E.2 Construction Code No. 1 

This construction code is designed mainly for thin-walled bars 
where the wall thickness only is to be changed during the optimization 
process. The size of the element, namely the cross-sectional area A, 
is then related to the moments of inertia by 



I- = j-A 
z z 


(E.2.1) 


vwhere i- and j- are the unit moments of inertia. 

-^y ^ z 

If P is positive (tensile), the element is redesigned with respect 
to the allowable tensile stress only. The stress ratio redesign formula 
(B.2.3) then becomes 


AVA = P/P* , (E.2.2) 

where A* is the improved area and P* = a*A, a* being the tensile strength. 

For negative (cpmpressive) values of P, Johnson’s parabolic formula, 
[Sec. Cl, Ref. 8] shown in Fig, E.2,1, is used as the failure criterion. 
For slender bars, failure is governed by Euler buckling, in which case 
the stress ratio redesign formula is 


AVA = -P/P^^ , (E.2.3) 

where = CJ^^A is the Euler buckling load of the current design. For 
short bars, the allowable stress is given by Johnsons parabola; the 
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= failure stress 
L = length of element 

r = radius of gyration of the cross section 

a* = compressive strength (yi^ld stress or crippling stress) 

^ 2 2 

Ocr = Ctt E/(L/r) = Euler buckling stress 

E = Young's modulus 

C = end-fixity coefficient for buckling 

Figure E . 2 . 1 

Johnson's Parabolic Formula 
resulting stress ratio formula can be shown to be 

PP 

A*/A = ~ p rp* -p*/4') (E.2,4) 

c^ cr c ^ 

where P* = a*A and P = a A. The ratio A' /A is computed about both 
c c cr cr ^ 

principal axes of bending, and the larger of the values is chosen as the 
stress ratio of the element. 

The elements may be guarded against torsional buckling (applicable 
to open sections), or local buckling of the wall by choosing appropriate 
minimum size constraints. 


I 
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VI B. Material Properties Data 

A separate data deck, described below, is required for each 
material. The decks do not have to be in a numerical sequence of 
material numbers. 

Material Control Card (215, FIO.O) 


1 

6 

11 20 

\_J_ 

NTC 

1 WT j 


N = Material number. 

NTC = Number of temperatures for which the properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight-density of the material. 

Material Properties Cards (5F10.0) 

One card for each temperature for which the properties of this material 
are given. Cards must be in ascending order of temperatures. 


1 

11 

21 

31 

41 50 

l_L 

E 

a 

°t* 

a* 

c 


PMAT(5) 5» 


PMAT(l) = Temperature for which the properties are given (T) . 

PMAT(2) = Young's modulus of elasticity (E) . 

PMAT(3) = Coefficient of linear thermal expansion (a). 

PMAT(4) = Tensile strength (a*) . 

PMAT(5) = Compressive strength (a*). If blank, computer sets PMAT(5)=PMAT(4) . 

VI C. Geometric Properties Data (15, 5X, 3F10.0) 

One card for each geometric property; cards do not have to be in 
numerical order. 


1 

6 

11 

21 31 40 

N 

blank 

AREA 

I- 1 I- 
y ' z 




^PGE0(2) >■ 
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N = Geometric property number. 

AREA = Cross-sectional area ox reference section (see below) . If 
blank, computer sets AREA = 1.0. 

PGEO(l) = Moment of inertia of the reference section about y-axis (I-). 

If blank, computer sets PGEO(l) = 10.0 **6, thus eli- ^ 
minating Euler buckling as a design consideration. 

PGE0(2) = Moment of inertia of the reference section about z-axis (I^). 

If blank, computer sets PGE0(2) = 10. 0**6^ eliminating ^ 
Euler buckling. 

The Geometric Properties Data is used solely for the computation of 
the unit moments of inertia j- and j-. Therefore, the reference section 
does not have to coincide witK the initial design, but can represent any 
acceptable intermediate design. The unit moments of inertia could, of 
course, be read in directly by setting AREA = 1.0 (or blank), in which 
case PGEO would be interpreted as the unit inertias. 

VI D. Element Load Multipliers (4F10.0) 

Four cards as shown below: 



Card 1: x-gravity 
Card 2: y-gravity 
Card 3: z-gravity 
Card 4: thermal loading 


EMUL = Fractions of the basic element loads (gravity loading in the 
three global coordinate directions and thermal loading) which 
are to be included in the element loads A, B, C and D. 

Element Load Multipliers simply define the element loads A, B, C and D. 
Various multiples of these element loads can be added to each structural 
load condition by the use of the Structural Load Multipliers (see Sec. D.l). 
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VI E.. Element Data (615, 4F10.0, 15) 

One card for each element; cards must be arranged in an ascending 
order of element numbers. 


1 

6 

11 

16 

21 

26 

31 

41 

51 

61 

71 75 

lEL 

II 

1 JJ 

IMAT j 

IGEO 

1 IDV 

FRC 

REFT 

ELPYY 

ELPZZ 

INC 


lEL = Element number. 

II = Number of node i (see Fig. E.1.1). 

JJ = Number of node j . 

IMAT = Material number of element, 

IGEO = Geometric property number of element. 

IDV = Design variable number of element. 

FRC = The design variable fraction r). in eqn. (B.1.2). If blank, 
computer sets FRC = 1.0. 

REFT = Reference temperature see Nodal Point Data in Sec. D.l. 

ELPYY = The end-fixity coefficient C for buckling about y-axis 

see Fig. E.1.2. 

ELPZZ =: The end-fixity coefficient C for buckling about z-axis. 

INC = Node number increment in autciaatic generation of element 
data (see below). If blank, computer sets INC = 1. 

Automatic element data generation if there exists a series of elements 

lEL = m, m+1, m+2,..., which satisfies the following requirements: 

(a) The node numbers form the sequences 
II = i, i+INC, i+2-INC,. . . 

JJ = j, j+INC, i+2*INC,...; 

(b) The rest of the data is identical for all elements of the series, 
then only the last card of the series will be required. The element data 
for the other elements in the series will be generated automatically. 
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E.3 Construction Code No, 2 

This construction code assumes that all the cross-sectional dimen- 
sions are changed by the same proportion upon redesign. Consequently, 
(E.2.1) is replaced by 


2 2 
I- = i-A , I- = j-A^ 
y -'y z •'z 


CE;3.1) 


This scheme allows the user to choose the optimal proportions of the 
, cross section (usually determined by local buckling considerations, 
such as crippling stress or torsional buckling), and assures him that 
these proportions are not changed upon redesign. 

If P is tensile, (E.2.2) remains valid as the redesign formula. 

For compression, the quadratic size-inertia relationships result in the 
following modifications to (E.2.3) and (E.2.4), respectively: 

. 1/2 


A’ /A = (-P/P^r)' 


(E.3.2) 


AVA = P*/(4P^P - P/P* . 


(E.3. 3) 


The input data and the remainder of the redesign procedure are 
identical to those used in construction code No. 1. 
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F. BEAM ELEMENT 


F.l General Information 



Figure F.1.1 

Typical Beam Element Showing Directions 
of Positive End-Forces and Moments . 

The element coor^iinate axes y and z shown in Fig. F.1.1, are the 
principal axes of the cross section. The direction of the y-axis is 
defined by the plane of i-j-k, where k is a "third” nodal point that 
must be specified in addition to the nodes i and j . The node k may be 
a conventional nodal point of the structure, or a point used solely for 
defining the principal axes. In the latter case, the motion code on 
the node data card should be used to eliminate all degrees of freedom 
for node k. 

Cubic polynomials in x are used for the lateral displacements, 
whereas the axial displacement is assumed to be linear as in the bar 
element. The bending moments are thus confined to linear functions 
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of X. It follows that the finite element analysis of beam elements is 
exact (apart from numerical errors) for elements that carry concentrated 
nodal forces or moments only. 

The basic element loads are gravity loads in the three global co- 
ordinate directions, and fixed-end forces referred to the element coordinate 
directions. Neither thermal stresses nor temperature-dependent material 
properties are included in the current version of the program. 

Distributed lateral loads are specified in the form of fixed-end 
forces and moments. A fixed-end force (or moment) is simply the end 
force (or moment) acting in a clamped-clamped beam subjected to the 
specified lateral loading. An example on the use ot fixed-end forces 
is shown in Fig. 1.2. 
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Figure F.1.2 

Example on the Use of Fixed-End Forces. 
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The use of fixed-end forces results in the correct end- forces 
(and moments) in a beam element, but the distribution of the internal 
forces (and moments) is generally inaccurate between the ends. 

Any degree of freedom of an end -node (such as i or j) can be 
connected to the degrees-of-freedom of a "master” node m. If such a 
connection is present, the number of the master node should be used as 
the motion code for the "slave" node (see Nodal Point Data in Sec. D.l). 
Note that this scheme precludes using node number 1 for the master node. 

If the x-displacement of node i is specified as a slave of node m, 
the displacement is expressed in terms of the degrees-of-freedom of the 
master node as follows; 


u . = u + (z .-2 )0 - (y.^y )0 

XI xm ^ 1 m^ ym -^m*' zm 


(F.1.1) 


The corresponding formulas for the y and z-displacements can be obtained 
from (F.1.1) by cyclic pernutation of x,y,z. 

If the rotation about x-axis of node i is specified to be the slave 
degree-of-freedom, the computer sets 


0 . = 0 
XI xm 


(F.1.2) 


Analagous equations are used for slave rotations about the y and z-axis. 

Slave-master relationship is used to represent various rigid 
connections between two nodes. The edge-reinforced plate in Fig. F.1.3 
can be offered as an example. By specifying all the degrees-of-freedom 
of node j to be the slaves of node m, a rigid link is created between 
the two nodes. The link enforces displacement continuity between the 
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Figure F.1.3 

Example of Master and Slave Nodes. 


beam and plate element, but at the same time it accounts for the inter- 
element eccentricity e. 

It is important to note that the eccentricity e is kept constant 
throughout the redesign procedure, i.e,, no allowance is made for the 
changes in e caused by redesign. 

The printed output following each analysis consists of the end- 
forces and moments shown in Fig. F.1,1. 

The first card of the element data deck for each construction code 
used must be: 

VI A. Element Control Card (615) 

1 6 11 16 21 26 50 

MTYPE I NUME j NUMMAT | NUMGEO | KODE |nUMFX 
-S NPAR(6) ^ 

NPAR(l) = Code for element type (MTYPE). For beam elements use NPAR(l) = 
NPAR(2) = Number of elements with the specified construction code (NUME). 
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NPAR(3) = Number of different materials used for the specified con- 
struction code (NUMMAT) . 

NPAR(4) = Number of different geometric (cross-sectional) properties 
used for the specified construction code (NUMGEO) . 

NPAR(5) = Construction code (KODE) . 

NPAR(6) = Number of different fixed-end force sets used for the specified 
construction code (NUMFX) . 

The rest of the element data is dependent on the specified construction 

code, and is listed separately below. 
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F. 2 Construction Code No. 1 

The size A of the element is taken as the cross-sectional area, and 
the size-stiffness relationships are assumed to be linear: 

I; ■ ' Y’ h ' ^5'" ■ 

where I- is the torsional constant of the cross section, I- and I- 
X * y z 

are the principal moments of inertia, and j-', j-, j- are the size- 
independent unit values . 

This construction code is intended primarily for the use of thin- 
walled sections where the wall thickness only is to be varied during 
design. Equations (F.2.1) represent a good approximation for closed 
sections, provided that the wall thickness is not too large. For open 
sections, the expressions for I- and I- are still valid, but the 

3 

torsional constant has the form I- = j-A . It must be noted, however, 
that I- is very small in comparison to I- and 1- for thin-walled, 
open sections. Therefore, its contribution to the rigidity of the 
structure can be neglected in most cases, i.e. j" “ ^ niay be used. 

The construction code is also valid for plane bending of beams 
with rectangular cross sections where the width of the cross section 
is subject to design. 

The stress ratio redesign formula is based on the allowable 
stresses only; instability criteria are not used. Local instabilities 
could be guarded against by choosing adequate minimum size constraints. 
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The stresses are calculated at four points at each end section 
of the element, as indicated in Fig. F. 2.1 by points A, B, A' and B’. 
Locations of A and B on the cross section are chosen by the user; 
points A* and B* are located by the computer to be in "mirror image" 
position of A and B. 



z 


Circular 
KSEC = 3 

Figure F . 2 . 1 

Allowable Families of Cross Sections. 
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The input is limited to three families of cross sections: 

(a) Cross sections with at least one axis of symmetry. It is also 
advisable to have the centroid and the shear center coincide, 
since the program takes the center of twist to be the centroid 
of the cross section. The axis of symmetry must be specified 
as the y-axis. 

(b) Zee section with identical flanges. The y-axis must be speci- 
fied as the principal axis that passes through the centroids of 
the flanges as shown in Fig, F.2.1. 

(c) Circular cross section. The orientation of y and z axes is 
completely arbitrary; the computer will automatically rotate the 
axes such that the y-axis will coincide with direction of the 
resultant moment M (see Fig. F.2.1). 

The user must specify the section moduli at points A and B only; the 

computer will automat;ically generate the moduli for A' and B’ as shown 

below. 

(a) Symmetric sections: 


'Va' 

= -(Va> 






= -CZ-)b. 
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Zee sections: 
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The normal stresses due to bending and axial load are computed 
from the formulas 

ag = ? (M-/Z- - M-/Z-) . (F.2.3) 

Oa = ^ R-/A , CF.2.4) 

respectively, where the minus sign applies to node i, and the plus sign 
to node j . The shear stress due to torsion is obtained from 

I = M-/Z- . CF.2.5) 

The stress ratio redesign formula is based on the modified Von Mises 
yield criterion 

+ (T/T*)^ = 1 , (F.2.6) 

where a = and a*, T* are the allowable stresses. If a > 0, 

o* is ‘taken as the allowable tensile stress a*; if a < 0, a* represents 
the allowable compressive stress o*. It should be noted that (F.2.6) 
reduces to the original Von Mises criterion if a* = a* and t* = a*/v^. 

Because only the wall thickness changes during design, the 
section moduli have the form 


Z- = z-A , 

X X 


Z“ = z-A , Z- = z-A , 
y y z ^ 


(F.2.7) 


where z-, z-, and z- are the size-independent unit section moduli. 
Utilizing (F.2.7) and (F.2.3) to (F.2.6) in the failure criterion (F.2.6), 
we obtain for the stress ratio redesign formula the following equation: 


A’ 

A 


“ V^B.2 
C^;r— ) + 



(F.2.8) 


where a^, a^, and x are the stresses of the current design A. 
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VI B. Material Properties Data (15, 5X, 6F10.0) 

One card for each material used; the cards do not have to be in numerical 
sequence. 

1 6 11 21 31 41 51 61 70 


N 

blank 

WT 




^ rMAl J - " 


N = Material number. 

WT = Weight-density of the material. 

PMAT(l) = Young’s modulus of elasticity (E) . 

PMAT(2) = Poisson’s ratio (v) . 

PMAT(3j = Allowable tensile stress (c?*) . 

PMAT(4) = Allowable compressive stress (a*). If blank, computer sets 
PMAT(4) = PMAT(3). 

PMAT(5) = Allowable shear stress (x*) . If blank, computer sets 
PMAT(5J = 0.577*PMAT(3) . 

VI C. Geometric Properties Data C2I5, 4F10.0/6F10.0) 

Two cards for each geometric property; cards do not have to be in 


numerical sequence. 


1 6 11 21 31 41 50 


N 

KSEC 

AREA 

I- I I- 

X 1 y 

I- 

Z 

1 

PGE0(9) H 


1 

11 

21 

31 

41 

51 60 

A 

.,A 

^A 




Z“ 

Z- 

Z- 

Z- 

Z- 

Z- 

X ' 

y 

z 

X 

y 

z 


-6 PGE0(9) H 


Card 1 


Card 2 


N 


= Geometric property number. 
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KSEC 

AREA 

PGEO(l) 

PGEO(2) 

PGEO(3) 

PGEO(4) 

PGEO(5) 

PGEOC6) 

PGEO(7) 

PGEOC8) 

PGEO(9) 


= Code for cross section (see Fig. F.2.1). 

KSEC = 1 or blank: s/minetric section. 

KSEC = 2: Zee section. 

KSEC = 3: circular section. 

= Cross-sectional area of reference section. 

If blank, computer sets AREA = 1.0. 

= Torsional constant of the reference section (I-). 

= Moment of inertia of the reference section about y-axis C^’) • 

= Moment of inertia of the reference section about z-axis (I-). 
If KSEC = 3 (circular section), computer sets ^ 

PGE0(3) = PGE0(2). 


= Section modulus of point A of the reference section for 
torsion (Z*-) . 

= Section modulus of point A of the reference section for 
bending about y-axis (Z^) . 

= Section modulus of point A of the reference section for 
bending about z-axis (Z^) . 

= Section modulus of point B of the reference section for 
torsion (z|) . 


= Section modulus of point B of the reference section for 
bending about y-axis (Z^) , 

= Section modulus of point B of the reference section for 
bending about z-axis (z|) . 


Note: If PGEO(i), i = 4, 5,..., 9 are left blank, the 

corresponding torsional or bending stress is set to zero 
in the analysis. For KSEC = 3 (circular section), stresses 
are calculated at points A and A* only; consequently 
PGE0(6) to PGE0(9) are set to zero by the computer. 
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The Geometric Properties Data is used only for the computation of unit 
inertias and section moduli* The reference section thus does not have 
to be the initial design. 


VI D. Element Load Multipliers (4F10.0) 
Three cards as' shown below. 


1 

11 

21 

31 40 

A 

B 

c 

D 

A 

B 

C 

D 

B 

B 

C 


ENIUL (3.4) 


Card 

1: 

x-gravity 

Card 

2: 

y- gravity 

Card 

3: 

z-gravity 


ENIUL = Fractions of basic element loads (gravity loads in the global 
coordinate directions) which are to be included in the element 
loads A, B, C and D. 

Lumped nodal forces are used to represent gravity loads, as was done for 
Bar Elements (see Sec. E.l). Fixed-end forces are not computed. 

The element load multipliers define the contribution of gravity to element 
loads A, B, C and D. Additional contribution to the element loads is 
made by the fixed-end forces; this information is specified separately by 
Fixed-End Force Data and Element Data. 

Specified multiples of element loads A, B, C and D can be added to each 
structural load condition by the use of Structural Load Multipliers 
(see Sec . D. 1) . 

VI E. Fixed-End Force Data (15, 6F10.0/5X, 6F10.0) 

Two cards for each set of fixed-end forces used in the analysis; cards 
do not have to be in numerical order. If no fixed-end forces are present 
(NUMFX = 0), omit all cards. 


1 6 16 26 56 46 56 65 


N 

R- 

R- 

R- 

M- 

M- 

M- 

X 

y_ 

7 

X 

y 

z 

blank 

R- 

R- 

R- 

M- 

M- 

M- 

X 

V 

z 

X 

J 

z 


-e SFT(12) H 


Card 1 : node i 

Card 2 : node j 
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N 

= Fixed -end force set number. 


SFT 

= Fixed-end forces and moments, 
in Fig. F.1.1. 

. Use the sign convention shown 

VI F 

Element Data (715, FIO.O, 415, 

216, 13) 


One card for each element; cards must be in an ascending order of 
element numbers. 


1 6 11 16 21 26 31 36 46 51 56 61 66 72 78 80 


lEL 

I 1 J i K 

IMAT 

IGEO 

IDV 

FRC 

a1 b1 cl D 

NI 1 NJ 

INC 

1 

IEC3)— ^ 


h^LC(4)-^ 

-^JC(I2) — H 



lEL = Element number. 

I = Number of node i (see Fig. F.1.1). 

J = Number of node j . 

K = Number of node k. This node defines the direction of y-axis, 

and should not lie on the line of i- j , 

IMAT = Material number of element. 


IGEO = Geometric property number of element. 
IDV = Design variable number of element. 


FRC = The design variable fraction m in eqn. (B.1.2). If blank, 
computer sets FRC = 1.0. 

LC(4) = Numbers of fixed-end force sets that are to be included in 

element loads A, B, C and D, respectively (also see Element 
Load Multipliers) . 

NI = End release code for node i. The end release code consists of 
a six digit number, each digit corresponding to an end-force 
as shown below. 

If any of the end-forces is known to 
be zero, due to the presence of a 
hinge or a roller, the digit one must 
be used; otherwise the digit should be 
left blank. 


1 2 3 4 5 6 


R- 

R- 1 

R- 

M- 1 

M- 

M- 

X 

_ZJ 

z 

x 

y 

z 
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NJ = End release code for node j . 

INC = Node number increment in automatic generation of element data 

(see below'). If blank, computer sets INC = 1. 

Automatic element data generation if there exists a series of elements 

lEL = m+1, m+2,..., which satisfies the following requirements: 

(a) The numbers of the end-nodes form the sequences 

l=i, i+INC, i+2*INC,..., 

J=j, j+INC, j+2*INC,...; 

(b) The rest of the data is identical to all elements of the series, 
including the node number K, 

then only the last card of the series will be required. The element data 
for the other elements of the series will be generated automatically. 
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F . 3 Construction Code No. 2 

In this construction code all the dimensions of the cross section 
are scaled by the same factor upon redesign. This scheme gives the 
user an opportunity to choose optimal proportions for the cross section, 
and maintain these proportions throughout the design process. 

Tlie expressions for the bending and torsional properties of the 
cross section now take form 


I- 


X 





Z- 


X 




Z- = 
z 



and the formula for stress ratio redesign becomes 


(F.3.1) 

(F.3.2) 



(F.3.3) 


This equation is solved iteratively for A* /A in the manner described 
below. 

Let r"^ = (AVA)*^^be the ratio obtained from the previous iteration 
(iteration v) , and r^^^ = (A’/A)^^^ the improved value predicted by the 
current iteration. If 


^ I r \2 

2 —It** " 


(F.3.4) 


the improved value is obtained from the solution of the cubic 
, v+1,1/2 

equation in (r ) : 




(F.3.5) 
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If (F.3.4) is not satisfied, the iteration formula is taken as cubic 
v+1 

equation in r : 


, V+l-3 V+1 



+ 




i<y*r J 


0 


(F.3.6) 


If T = 0, (F.3.5) yields an exact solution of (F.3.3). 

Similarly, (F.3.6) is an exact redesign equation if = 0 

or a„ = 0. In both cases only one iterative cycle is needed. Other- 

D 

wise the number of iterations is limited to six in each redesign cycle, 
i.e., the results of the sixth iteration are used for the new design 
even if the iterative solution failed to reach convergence. 

The input data and the remaining details of analysis and re- 
design procedures are identical to those used for Construction Code 


No. 1. 


G. PLANE STRESS QUADRILATERAL OR TRIANGULAR ELEMENT 


G.I General Information 



(a) 



(b) 


Figure G.1.1 


Quadrilateral and Triangular Plane Stress Elements. 



I 
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The basic plane stress element is the quadrilateral shown in 
Fig. G.l.l(a). The nodes 1 and 2 determine the local Cartesian coordinate 
axes X and y, and the relative orientation o£ the natural element 
coordinates s and t in the manner indicated on the drawing. Note that 
the natural element coordinates take the values +1 or -1 at the sides of 
the element. All four nodes of the element should be on the same plane. 

If the nodes are not co-planar errors will occur in the computation of 
element properties. The magnitude of the error becomes more serious 
with increasing warp of the element. 

The displacement field within the element is taken in the form 

4 

u(s,t) = I h^(s,t)u^ + h^(s,t)a^ + h^(s,t)a 2 , 
i=l 

4 

v(s,t) = Z h^(s,t)v^ + hg(s,t)a 2 + h,(s,t)a. , (G.1.1) 

i=l 

4 

w(s,t) = Z h. (s,t)w , 
i=l 

where u, v and w are the components of the displacement field in the 
direction of the local Cartesian coordinate axes, u. , v. and w. denote 
the displacement components of node i, and are generalized coordinates 
in addition to the nodal displacements. The shape functions in (G.1.1) 
are 

hj = i(l-s)(l-t) = i(l+s)(l-t) 

h 3 = i-(l+s) (1+t) = i(l-s) (1+t) (G.1.2) 

hs = i(l-s)2 hg = i(l-t)2 . 
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The functions to are Lagrangian interpolation polynomials; they 
result in a linear displacement field along each side of the element, 
and thus maintain displacement continuity between the elements. The 
functions h^ and h^, on the other hand, represent incompatible displace- 
ment modes, because they cause the sides of the elements to be deformed 
into parabolas . 

According to the examples in [2] , the introduction of the incompatible 
modes can lead to a significant improvement of the prebuckling analysis 
for certain problems. 

Note that h^ and h^ vanish at the nodes. Consequently, ctj’s are 
internal degrees of freedom, which can be eliminated at the element level 
by the static condensation procedure. The user of the program, however, 
has the option of suppressing the incompatible modes altogether, in which 
case the computer would set = 0, j = 1,2, 3, 4. 

All the area integrals in the computation of element properties are 
carried out numerically, using Legendre-Gauss quadrature with four 
integration points. It can be shown that the numerical integral is 
exact only for rectangular elements. Since the integration error in- 
creases with the skewness of the element, highly skewed quadrilaterals 
are to be avoided if possible. 

The triangular element used in the program is simply a limiting 
case of the quadrilateral when two of the nodes are allowed to coincide, 
as seen in Fig. G. 1.1(b). It requires no special provisions within 
the program. 

Tlie basic element loads are gravitational forces in the three 
coordinate directions, thermal loading, and in-plane pressure acting on 
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one side of the element (the pressure loading is not allowed for the 
undirectionally reinforced element described under Construction Code 1). 

Temperature-dependent material properties may be used in the pro- 
gram, in which case the material properties should be specified for several 
temperatures. The material properties of each element are obtained 
from the specified properties by linear interpolation. The highest and 
lowest temperatures at which the material properties are prescribed should 
cover the entire range of element temperatures. 

The output of each analysis cycle consists of the in-plane stress 

resultants at the origin of the natural coordinates point 0 in 

Fig. G.l.l(a), and at points A, B, C, and D. The stress resultants at 
point 0 are printed out twice: once with respect to the local Cartesian 

coordinates x and y, and once with respect to the coordinates x and y 
defined by the angle 6 in Fig. G.1.1. The mid-side forces are expressed 
in terms of the directions x^, y^, etc. shown in Fig.. G.1.1 (a). The 
convention for positive stress resultants is given in Fig. G.1.2. 



Figure G.1.2 


Positive Stress Resultants 
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By using an appropriate value of the stress-printout option code NS, 
the user can suppress the computation of some or all of the stresses. 

The element data deck for each construction code must start with: 

VI A. Element Control Card 



NPAR(l) = Code for element type (MTYPE) . For plane stress elements use 
NPAR(l) = 3. 

NPARC2) = Number of elements with the specified construction code (NUME) . 

NPAR(3) = Number of different materials used for the given construction 
code (NUMMAT) . 

NPARC4) = Maximiom number of temperatures for which material properties 
are given (NUMTC) . If blank, computer sets NPAR(4) = 1, 

NPAR(5) = Construction code (KODE) . 

NPAR(6) = Code for use of incompatible displacement modes (n) . 

NPAR(6) = 0: use incompatible modes. 

NPAR(6) > 0: suppress incompatible modes. 

NPAR[7) = Number of genetic properties (not used for Construction 
Code No. 2) . 
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G.2 Construction Code No. 1 

This construction code deals with unidirectional ly stiffened panels. 
All the cross-sectional dimensions shown in Fig. G.2.1, including the 
stiffener spacing w, are assumed to change by the same proportion during 
redesign. 

w 



Note: A and I are the cross-sectional area and moment of inertia of 

s s 

the stiffener, respectively. 

Figure G.2.1 

Cross Section of Unidirectionally Stiffened Panel. 


The sheet thickness t is chosen as the size of the element, i.e., it 
is taken as the independent dimension of the cross section. 

For the purposes of analysis, the panel is treated as a homogeneous, 
orthotropic slab with equivalent extensional and shear flexibilities: 



1/t 

-v/i 


-v/i 

1/t 



% 


y 


xy 


(G.2,1) 


0 


0 2(l+v)/t 
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In (G.2.1), X and f are the principal material coordinates shown in 

Fig. G.2.2; N^, and N/v/x represent the average membrane stress 
X y xy 

resultants; and 

t = t + A^/w 

is the "average” thickness of the panel. 



Figure G.2.2 

Unidirectionally Stiffened Panel 


The following possible modes of failure are considered in redesign 
operation: 

(1) stresses exceeding their allowable values in the sheet; 

(2) general buckling of the panel; 

(3) buckling of the sheet between the stiffeners; 

(4) failure of stiffeners. 

The stress resultants acting at the middle of the panel (s = 0, t = 0) 
only are used in evaluating these failure modes. 
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Other types of local buckling failure can be handled simply by 
choosing suitable cross-sectional proportions of the panel, i.e., pro- 
portioning the reinforcement such that sheet buckling would always occur 
prior to (or simultaneously with) the local buckling. The user may even 
optimize the cross-sectional proportions by maximizing i:he configuration- 
efficiency coefficient of the panel in a manner similar to Ref. [9], p. 107. 

(1) Stress limits 

The modified Von Mises yield criterion 

= 1 (G.2,2) 

is used for the stress-constrained design of the sheet, where a* and t * 
are allowable normal and shear stresses, respectively. Different allowable 
stresses may be specified for tension (a*) and compression (a*) . Equa- 
tion (G. 2. 2) reduces to the original Von Mises yield criterion only if 

a* = a* and t* = a*/y/Z, 
t c 

Using (G.2.2), the stress ratio redesign formula can be shown to be 

t'/t = [(N^/Np^ + CN./Np^- N^N./(N|Np + (N^./N|.)^] (G.2.3) 

where - a*t, NS. = a*t and NS^ = T*t . 

X " y xy 

(2) General buckling of panel 

The panel is treated as a homogeneous, orthotropic plate with simple 
supports. It is assumed that only the compressive stress in the 
direction of the stiffeners influences buckling. 

Since the buckling formula for a general quadrilateral is not 
available, the panel is treated as a rectangle with -the ’’effective” width 
b shown in Fig. G.2.2. The effective width is computed automatically; 
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however, the user has the option o£ substituting his own value. 

According to Ref. [8], Sec. C.2, p. 26, the compressive buckling 
stress of a stiffened plate is 


^ X cr 


kTT^E ^t^2 
2, '■b'' 


12(l-v‘-) 


(G.2.4) 


where the buckling coefficient is given by 


1 + 
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N wt 


(G.2.5) 


In (G.2.5), N is the number of bays between stiffeners , and 


2 

12 (iV) 

The remaining symbols are defined in Fig. G.2.1. 

For large numbers of stiffeners , we can approximate (N-l)/N = 1 and 
N = b/w, and obtain for the critical stress resultant 


where 


k, = 


(¥cr 


2(^)^ 


kjir^E 


2 4 

12(l-v )b 


(G.2.6) 




12(l-v")I 


wt' 


1 + 


A d^/I 
s s 


0.12+0.88t/t/ J 


1/2 


+ 1 
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Note that the coefficient k^ does not change upon redesign since the 
cross-sectional proportions are kept constant. The stress ratio re- 
design formula hence becomes 


tvt . 


(G.2.8) 
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(3) Buckling of sheet between stiffener s 



Figure G.2.3 

Sheet Between Stiffeners. 


The sheet between stiffeners , shown in Fig. G.2.3, is treated as a 
long, simply supported plate of effective width w^. The general form 
of the compressive buckling stress for each of the three load conditions 
shown is 


a = k 
cr 


TT E 

' t ' 

12(l-v^) 

•W 

i e/ 


(G.2.9) 


For simple supports, the buckling coefficients for longitudinal 
compression, transverse compression and shear are, respectively 

k^ = 4.00, k^ = 1.00, k^^ = 5.35. (G.2.10) 

X y xy 

If all loads are applied simultaneously, the buckling criterion 
(interaction formula) is taken as 

9 2 

-C^/COxs) - O/s/CO/s) + 

X ^ x*^cr y y"^cr xy xy"^cr 


1 


(G.2.11) 
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The interaction fomiula has an exact theoretical basis only when = 0 or 

y 

= 0; otherwise is an approximation. 

The stress ratio redesign formula obtainable from (G.2.9) to (G.2.11) 
is 


- N. - 4(t/t)N^ + {[N. + 4(t/t)N^]^ + [1.495(t/t)N..]^}^/^ 


2(N.) 

X cr 


, (G.: 


where 


(N^) = t . 

^ x'^cr x^cr 


Equation (G.2.12) remains valid if and N/s are positive (tensile) 

X y 


(4) Failure of stiffeners 



Figure G.2.4 

Stiffener and Sheet Treated as a Column. 


A stiffener and the attached sheet (see Fig. G.2.4) are treated as a 
column of length The value of ”a” is calculated by the computer 

in the manner indicated in Fig. G.2.2, or it may be user-specified. 


. 12 ) 
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If is positive (tensile), the failure analysis is skipped. For 
compressive loading, the Euler-Johnson failure criterion shown in 
Fig. E.2.1 is used. 

The cross-sectional area of the column in Fig. G.2.4 is 

A = + wt (G.2.13) 

and its moment of inertia can be shown to be 

I = wt^/12 + Ig + d^[wt(Ag/A)^ + A^(l-Ag/A)^] . (G.2.14) 

The unit moment of inertia, which is invariant of the value of t, 
can be calculated from 

j = I./t'* . (G.2.15) 


For long columns that fail by Euler buckling, the critical compressive 
stress resultant is 


_ Ctt^EI _ CTT^Ej ^3 

cr ■ 2 2. . , ’ 

aw a (w/tj 

from which we obtain the redesign formula 


(G.2.16) 


t'/t 


(-N^/N ) 

^ X cr 


1/3 


(G.2.17) 


Johnson's parabola, which governs the failure of short columns, 
leads to the following quadratic equation in t'/t: 

N*(t'/t)^ + N.(t/t') - (N*)2/(4N^P = 0 , (G.2.18) 

where N* = cr*t, a* being the compressive strength of the stiffener 
(yield stress or crippling stress) . 



G.2.8 


VI B Material Properties Data 

A separate data deck is required for each material used. The 
decks do not have to be in numerical sequence of material numbers. 
Temperature - Indep endent Pat a (215, FIO.O) 


1 6 11 20 


NTC 

WT 


N = Material number 

NTC = Number of temperatures for which properties of this material 
are given. If blank, computer sets NTC = 1. 

IVT = Weight-density of the material. 

Temperature-Dependent Data (8F10.0) 

One card for each temperature for which the properties of this 

material are given. Cards must be in ascending order of temperatures. 


1 

11 

21 

31 

41 

51 

61 71 80 

[I 

1 E 

1 V 

i a 

1 a* 

‘ t 

1 a* 
* c 

1 T* 1 a* 
‘ ' s 


PMAT(8) ^ 


PMAT(l) 

PMAT(2) 

PMAT(3) 

PMAT(4) 

PMAT(5] 

PMAT(6) 

PM^\T(7) 

PMAT(8) 


Temperature for which the properties are given (T) . 

Young *s modulus of elasticity (E) . 

Poisson’s ratio (v) . 

Coefficient of linear thermal expansion (a) . 

Tensile strength of sheet C*^*) • 

Compressive strength of sheet (o'*). If blank, computer sets 

^ > PMAT(6) = PMAT(5), 

Shear strength of sheet (t*). ^ PMAT(7) = 0.577 PJ4AT(5). 

Compressive strength (yield or crippling stress) of stiffener 
(a*). If blank, computer sets PMAT(8) = PMAT(6). 
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VI C Geometric Properties Data (15, 6F10.0) 

A separate card is required for each element geometry. The cards 
do not have to be in numerical sequence of geometric property numbers. 


1 

6 

16 

26 

36 

46 

56 65 

[n 

TH 

w 

SA 

SI 

D 

WE 


= Geometric property number. 

= Sheet thickness t of the reference cross section (see note below) 

= Distance w between the stiffeners of the reference cross section. 

= Cross-sectional area A of the stiffener in the reference cross 
section, 

= Moment of inertia Ig of the stiffener about its own centroidal 
axis in the reference cross section. 

= Distance d between the centroid of the stiffener and the mid-plane 
of the sheet in the reference cross section. 

= Effective width Wg of the sheet between the stiffeners in the 

reference cross section (used in buckling analysis of the sheet). 
If blank, computer sets WE = W. 

The Geometric Properties Data is used only for the computation of 
unit cross-sectional properties, i.e., properties for unit sheet thickness 
The reference cross section, therefore, can be chosen as any design 
that has the desired cross-sectional proportions; it does not have to 
be the initial design. 

VI D Element Load Multipliers (4F10.0) 

Four cards as indicated below: 


N 

TH 

W 

SA 

SI 

D 

WE 


1 

11 

21 

31 40 

thermal 

x-gravity 

y-gravity 

;:-gravity 

thermal 

x-gravity 

y-gravity 

z- gravity 

thermal 

x-gravity 

y-gravity 

z- gravity 

thermal 

x-gravity 

y-gravity 

2 -gravity 




Card 1: Element load A 
Card 2: Element load B 
Card 3: Element load C 
Card 4: Element load D 
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E^^UL = Fractions of basic element loads (thermal loading and gravity 
loads in the three global coordinate directions) which are to 
be included in the element loads A,B,C and D* 

Element Load Multipliers simply define the element loads A,B,C and D. 

Various multiples of these element loads can be added to each structural 

load condition by the use of Structural Load Multipliers (see Sec. D,l). 

VI E Element Data (715, 5X, 4F10.0/2F10.0, 215) 

Two cards for each element; cards must be arranged in an ascending 

order of element numbers. 


1 6 11 16 21 26 31 36 41 51 61 71 80 


lEL 

I i J 1 K 

L 

IMAT 

IDV 

blank 

FRC 

REFT 

AA 

AB 


^ IE(4) H 



1 11 21 26 30 


BETA 

EFC 

NS 

INC 


Card 1 


Card 2 


lEL = Element number, 

IE = Node numbers I,J,K and L of the four corner nodes (see Fig. G.1.1). 

Nodes 1 and 2 should define the side tliat is most parallel 
to the stiffeners as indicated in Fig. G.2.2. For a triangular 
element (their use is not recommended for stiffened panels) , 
the same node number must be used for nodes 3 and 4. 

IMAT = Material number of element. 

IDV = Design variable number of element. 

FRC = The design variable fraction r). in eqn. (B.1.2): A. = p.D . 

Note that A. , the size of the element, is taken as the ^ ^ 
sheet thickness t. If blank, computer sets FRC = 1.0. 

REFT = Reference temperature (temperature of the stress-free state) 
of the element. 

AA = The effective panel length ”a” used in computing the Euler 

buckling load (G.2.16) of the stiffeners. If blank, computer 
calculates ^^a'‘ in the manner shown in Fig. G.2.2. 

AB = The effective panel width b used in computing the general 

buckling stress (G.2.6) of the panel. If blank, computer calcu- 
lates b in the manner shown in Fig. G.2.2. 
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BETA = The angle 6 between the stiffeners and the local x-axis 
(see Fig. G.2.2) . 

EFC = The end-fixity coefficient C used in Euler buckling load 

(G.2.16) of the stiffeners. If blank, computer sets EFC = 1.0 
(simple supports) . 

NS = Stress printout code: 


NS 

Points for which stresses are printed 
(see Fig. G.1.1) 

1 

None 

3 

0 

6 

0,A 

12 

0,A,B,C 

15 

0,A,B,C,D 


If blank, computer sets NS = 3. If NS = 15 and the element is 
triangular, computer sets NS = 12. Note that design of the 
element is always based on stresses at point 0, regardless of 
which printout code is used. 

INC = Node number increment in automatic generation of element data 
(see below). If blank, computer sets INC = 1 . 

Automatic element data generation if there exists a series of ele- 

ments lEL = m, m+1, m+2,..., which satisfies the following requirements: 

(a) the node numbers form the sequences 
I=i, i+INC, i+2*INC,. . 

J=j, j+INC, j+2*INC,..., 

K=k, k+INC, k+2*INC,. . 
il+INC, £+2*INC,..., 

(b) the remainder of the data is identical for all elements of the 
series , 

then only the last card of the series will be required. The element 
data for the other elements of the series will be generated automatically. 



G.3.1 


G.3 Construction Code No. 2 

Isotropic, homogeneous panels subjected to stress constraints 
only are treated under this construction code. The variable dimension 
of the element (element size) is the panel thickness t. 

Two failure criteria are used in the stress ratio redesign: 

a) Von Mises yield criterion: 

- (.a^/ o*) o*) = 1 , (G.3.1) 

where and are the principal stresses, and a* is the allowable 
normal stress. Different allowable stresses may be specified for tension 
(a*) and compression (a*) . The corresponding stress ratio redesign 
formula is 

t’/t = [(Nj/N*)^ + (NyN*)^ - (N^/N*)(N2/N*)]^/^ , (G.3.2) 

where N* = a*t, and N^, are the principal stress resultants. 

b) Maximum shear stress theory of failure: 


V/T*=l , (G.3.3) 

T being the maximum shear stress and t* its allowable value. This 
max ^ 

leads to the stress ratio redesign formula 


tVt = (N ) /N* 

' s^max s 


(G.3.4) 


where (N ) is the maximum shear stress resultant, and N* = x*t. 

^ s^max s 

The value of t*/t is chosen as the larger of (G.3.2) and (G.3.4); 
however, the user can suppress the use of (G.3.4) by setting T* = 0 
on the material data cards. 
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The redesign formulas are applied to all points of the panel where 
the stress printout is requested (points 0,A,B,C or D shown in Fig. G.1.1). 
This is in contrast to Construction Code No. 1, where point 0 only was 
used for redesign. 

VI B Material Properties Data 

A separate data deck is required for each material used. The decks 
do not have to be in numerical sequence of material numbers. 

Material Control Card (215, FIO.O) 


1 

6 

11 21 

N 

NTC 

1 WT 


N = Material number. 

NTC = Number of temperatures for which properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight-density of the material. 

Temperature-Dependent Data (7F10.0) 

One card for each temperature for which the properties of this 

material are given. Cards must be in ascending order of temperatures. 


1 11 21 31 41 51 61 70 


iim 

E 

V 

a 

mm 

o* 

im 





mmM 

c 



H PMAT(7) H 

PMAT(1)= Temperature for which the properties are given (T) . 
PMAT(2)= Young’s modulus of elasticity (E) . 

PMAT(3)= Poisson’s ratio (v) . 

PMAT(4)= Coefficient of linear expansion (a). 

PMAT(5)= Tensile strength (cr^ ♦ 
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PMAT(6)= Compressive strength (a*). If blank, computer sets PMAT(6)=PMAT(5) • 

PMAT(7)= Shear strength (t*) . If blank, maximum shear stress theory of 
failure will not be used in redesign. 

VI C Element Load Multipliers (4F10.0) 


Four cards as shown below: 


1 

11 

21 

31 

41 50 

thermal 

pressure 

x-gravity 

y-gravity 

z-gravity 

thermal 

pressure 

x-gravity 

y-gravity 

z -gravity 

thermal 

pressure 

x-gravity 

y-gravity 

z-gravity 

thermal 

pressure 

x-gravity 

y-gravity 

z-gravity 

-iS EMUL(4,5) 


Card 1 : A 

Card 2 : B 

Card 3 : C 

Card 4 : D 


EMUL = Fractions of basic element loads (thermal loading, pressure 
acting on side of element, and gravity loads in the three 
global coordinate directions) which are to be included in the 
element loads A,B,C and D. 

Element Load Multipliers define element loads A,B,C and D. Various 
multiples of tliose loads can be added to each structural load condition 
by the use of Structural Load Multipliers (see Sec. D.l). 

VI D Element Data (715, F5.0, 3F10.0, 215) 

One card for each element; cards must be arranged in ascending 
order of element numbers. 


1 

6 11 16 21 

26 

31 

36 

41 

51 

61 

71 

76 80 

lEL 

I 1 J 1 K 1 L 

IMAT 

IDV 



FRC 

REFT 

PRESS 

BETA 

NS 

INC 


IE(4) H 



lEL 


= Element number. 
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IE 

IMAT 

IDV 

FRC 

REFT 

PRESS 

BETA 

NS 

INC 


= Node numbers I,J,K and L of the four corner nodes (see 

Fig. G.1.1). Nodes 1 and 2 should define the side on which 
pressure load is applied, if any. For a triangular element, 
the same node number must be used for nodes 3 and 4. 

= Material number of element. 

= Design variable number of the element, 

= The design variable fraction in eqn. (B.1.2): A^ = 

Note that Aj^, the size of the element, is taken as the panel 
thickness t. If blank, computer sets FRC = 1.0. 

= Reference temperature (temperature of the stress-free state) 
of the element. 

= Compressive force per unit length (pressure resultant) applied 
to side 1-2 of the element. 

= The angle P which defines the x and y axes as shown in Fig. G.1.1. 

The angle is used only to control the stress printout at point 0 
(see Sec. G. 1) . 

= Stress printout code (see Sec. G.2). Note that the design of 
the element is based on the stresses at all the points for which 
stress printout is requested. 

= Node increment in automatic generation of element data (see Sec. G.2). 
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H, QUADRILATERAL SHEAR PANEL 



A shear panel, shown in Fig. H.1.1, is assumed to resist only shear 
tractions applied to its edges; it has no rigidity with respect to normal 
edge tractions. The resultant forces of the shear tractions Q^, i = 1,2, 3, 4, 
must form a self-equilibrating system, so that only one of resultants 
is statically independent. 

All four corner nodes should be on the same plane, otherwise errors 
will arise in the computation of element properties. 

Following NASTRAN [5] , the equivalent nodal forces F^ and F^ are 
obtained by "lumping” one half of each adjacent edge force into the corner. 
This method was chosen purely on the basis of convenience, because it 
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results in comer forces that are col linear with the diagonals of the 
element. Again, only one of the comer forces is statically independent, 
the remainder being determined by equilibrium requirements. 

In order to calculate the strain energy (i.e., the stiffness matrix) 
of the element, an assumption must be made regarding the distribution of 
the shear stress r in the panel. We adopted the distribution suggested 
by Garvey [10], which satisfies equilibrium equations, but not compati- 
bility of the resulting displacement field (except for a parallelogram 
and a rectangle). The assumed shear planes are shown in Fig. H.1.1; the 
magnitude of t is taken to be inversely proportional to the square of the 
distance d from the %aseline^' AB. The details of formulating the element 
properties, including the geometric stiffness matrix, can be found in 
Ref. [5]. 

Two potential difficulties arise in the use of shear panels: 

a) The state of stress assumed by Garvey can exist only in rectangu- 
lar panel, i.e., only a rectangle can be in a state of pure shear when 
subjected to tangential boundary tractions. It is important, therefore, 
not to use severely skewed shear panels in order to avoid erroneous 
results . 

b) The absence of extensional rigidity will, in general, lead to 

a singular stiffness matrix of the structure, unless each shear panel 
is surrounded by elements that are capable of absorbing the extensional 
comer forces (forces in direction other than and ^2^* elements 

are commonly used for this purpose, where the stiffness of each bar 
usually includes the extensional rigidities of the adjacent shear panels. 

Shear panels have traditionally been used in aerospace industry 
to idealize stiffened panels: the shear resistance of the sheet is 
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accounted for by the shear panel, whereas the extentional rigidities 
of the sheet and the stiffeners are lumped into surrounding bar elements. 

This idealization was undoubtedly useful in the days of hand compu- 
tation, but its utility has become dubious with the arrival of the finite 
element method. An orthotropic, plane stress quadrilateral element, such 
as described in Sec. G.2, is not only a more precise representation of 
a stiffened panel, but it also requires less input data and no more computer 
time. 

Shear panel has not been excluded from DESAP 2, partly because it 
may still be a useful element in special applications, and partly in 
recognition of its popularity with aircraft designers. 

Only gravity loads are pemiitted as the basic element loads. 

Thermal stesses are excluded, of course, since thermal expansion cannot 
be resisted by a state of pure shear.. Temperature-dependent material 
properties, however, are allowed. 

The printout that follows each analysis cycle consists of the 
shear flow each of the four nodes, and the average shear 

flow of the panel, given by 

1 4 

The first card of the element data deck for each construction code 
must be: 

VI A Element Control Card (515) 


1 

6 

11 

16 

21 25 

MTYPE 

AHTMC 

NU^^MAT 

NUMTC 

KODE 


NPAR (5) ^ 
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NPARCl) 

NPAR(2) 

NPAR(3) 

NPAR(4) 

NPAR(5) 


Code for the element type (MTYPE) . For shear panel use MTYPE = 4. 

Number of elements with the specified construction code (NUME) . 

Number of different materials used for this construction code 
(NUMMAT) . 

Maximum number of temperatures for which the material properties 
are given (NUMTC) . If blank, computer sets NPARC4) = 1. 

Construction code (KODE) . 
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H.2 Construction Code No. 1 

Homogeneous, isotropic shear panels are treated under this 
construction code. The thickness t of the panel is the variable dimension 
(size of the element) . Two design criteria are used in the redesign 
equations : 

(1) The average shear stress should not exceed its prescribed 
limit T*. The corresponding stress ratio redesign formula is 


tVt = N^yN* , (H.2.1) 

where is the average shear flow defined in (H.1.1), and N* = T*t. 

(2) The average shear stress should not exceed the buckling stress 
For a rectangular panel, the critical shear stress is 


T 

cr 


= k 




12(1-V ) 



(H.2. 2) 


where b is the unsupported width (shorter dimension) of the panel, and 
k represents the buckling coefficient for shear. The stress ratio re- 
design equation obtainable from (H.2. 2) is 

tVt = » (H.2.3) 


where N = T t. 
cr cr 

A good approximation for the buckling coefficient is 


k = + C2(b/a) , 0 < b/a < 1 , 


(H.2.4) 


where a is the longer dimension of the panel, and and are constants 
that depend on the edge support conditions only. The values of these 
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constants for various edge supports are given in Fig. H.2.1. The figure 
also shows the corresponding support codes ISU that are used to specify 
the boundary conditions on element cards. 


= 5.35 
= 3.99 
ISU = 1 


1 


V C = 5.35 
^ C = 7.25 
^ ISU = 4 


= 8.99 
=5.72 

I 2 


yy//;/y I SU = 2 


= 5.35 
= 5.63 
ISU = 5 


C, = 8.99 


77 "/T7 y/ T /rry/y y 


^1 

= 3.29 
ISU = 3 




7.07 

3.91 


ISU = 6 


Simply supported edge 


Clamped edge 


Figure H.2.1 

Edge Supports for Buckling of Shear Panels. 


Equation (H.2.4) agrees precisely with the theoretical buckling 
coefficients in the limiting cases b/a = 0 and b/a = 1; in between these 
values, the approximation is somewhat on the conservative side. 

The unsupported dimensions of the panel, a and b, are user-specified. 
However, if the dimensions are left blank, the computer will use the 
’’average” panel dimensions shown in Fig. H.2.2. 
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Figure H.2.2 

Shear Panel and Its Average Dimensions Used 
in Buckling Analysis. 


VI B Material Properties Data 

A separate data deck is required for each material used. The 
decks do not have to be in numerical sequence of material numbers. 
Material Control Card (215, FIO.O) 


1 

6 

11 20 

N 

NTC 

WT 


N = Material number. 

NTC = Number of temperatures for which the properties of this material 
are given. If blank, computer sets NTC = 1. 

IVT = Weight-density of the material. 

Temperature-Dependent Data (4F10 , 0) 

One card for each temperature for which the material properties are 

given. Cards must be in ascending order of temperatures. 
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1 

11 21 31-40 

T 

E 1 V 1 T* 

PMAT(4) H 


PMAT(l) = Temperature for which the material properties are given (T) • 
PMAT(2) = Young’s modulus of elasticity (E) • 

PMAT(3) = Poisson’s ratio (v) . 

PMAT(4) = Shear strength (T*) . 

VI C Element Load Multipliers 
Three cards as shown below. 


1 

11 

21 

31 40 

m 

B 

C 


A 

B 

C 

D 

m 

B 

c 

D 


H EMUL(4,4) 


Card 1 : x-gravity 
Card 2: y-gravity 
Card 3; z-gravity 


EMUL = Fractions of basic element loads (gravity forces in the three 
global coordinate directions) which are to be included in 
the element load cases A,B,C and D. 

Element Load Multipliers define element loads A,B,C and D, Various 
multiples of these loads can be added to each structural load condition 
by the use of Structural Load Multipliers (see Sec. D.l). 

VI D Element Data ( 15, 3F10.0, IS) 

One card for each element; card must be in ascending order of ele- 
ment numbers . 


1 

6 11 16 

21 

26 

31 

36 

41 

51 

61 

71 75 

[ lEL 

I 1 J 1 K 

1 L 

IMAT 

IDV 

ISU 

FRC 

AL 

BL 

INC 

i 

K— IE(4) — 

>• 



lEL = Element number. 
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I 

J 

K 

L 

IMAT 

IDV 


ISU 


FRC 


AL 


BL 


INC 


= Node number of node 
= Node number of node 
= Node number of node 
= Node number of node 


1 (see Fig. H. 2 , 2) , 

2 . 

3. 

4. 


= Material number of element. If blank, computer sets IMAT = 1. 

= Design variable number of element; 

= Code for edge support conditions used in lodal buckling analysis 
(see Fig. H.2.1). If blank, local buckling analysis will be 
suppressed. 

= The design variable fraction in eqn. (B.1.2): Aj^ = 

Note that Aj^, the size of the element, is taken as the panel 
thickness t. If blank, computer sets FRC = 1.0. 

= The longer unsupported dimension "a" of the panel used in local 
buckling analysis. If blank, *'a** is calculated as shown in 
Fig. H.2.2. 

= The shorter unsupported dimension ”b*' of the panel used in local 
buckling analysis. If blank, **b” is calculated as shown in 
Fig. H.2.2. 

= Node number increment in automatic generation of element 
data (see below). If blank, computer sets INC = l. 


Automatic element data generation if there exists a series of 


elements lEL = m, m+1, m+2,... which satisfies the following require- 
ments ; 

a) the node numbers form the sequences 

I = i, i + INC, i+2*INC,. . . , 

J = j , j+INC, j+2*INC,. . . , 

K = k, k+INC, k+2*INC,... . . 

L = Z, £+INC, £+2*INC,. . . ; 

b) the remainder of the data is identical for all elements of the 
series , 

then only the last card of the series will be required. The element 


data for the other elements of the series will be generated automatically. 
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I. PLATE QUADRILATERAL* OR TRIANGULAR ELEMENT 
I.l General Information 



Figure I. i.l 


Quadrilateral Plate Element. 

The basic plate-shell element is the plane quadrilateral shown in 
Fig. 1.1.1. Tlie element is made up of four sub-triangles; the coordinates 
of the internal node (node 5) are obtained by averaging the coordinates 
of the four corner nodes.. The local Cartesian coordinates x and y are 
determined by the mid-points of the sides in the manner shown on the 
drawing. If the element is orthotropic, x and y denote the principal 
material coordinates; their direction is specified by the angle 3. 

The membrane behavior of the element is obtained by treating each 
sub-triangle as a constant strain element, and then eliminating the 
degrees -of -freedom of node 5 by matrix condensation. 

The procedure of Clough and Felippa[ll] is employed to derive the 
bending properties of the quadrilateral. A so-called LCCT-9 element, 
with nine degrees of freedom, is used for each sub triangle. The 
lateral displacement field of this element is a cubic function, but the 
displacements are partially constrained such that the normal slope along 
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each side is linear. After the four triangles are assembled into the 
quadrilateral element, the degrees of freedom of node 5 are again 
eliminated by static condensation. 

In calculating the contribution of the lateral displacements to 
the geometric stiffness matrix of the element, one modification is made 
to the above procedure: the shape functions used by Clough and Felippa 

for the subtriangles are replaced by the nonconforming functions 

suggested by Zienkiewicz Ref. [12], equation (10.26). A quint ic 

numerical integration formula (p. 151 of Ref. [12]) is then used to 
compute the unit geometric stiffness matrices for each of the subtriangles. 

The membrane and bending characteristics of the quadrilateral are 
condensed and stored separately in the form of unit stiffness matrices 
and load vectors as described in Sec. B.l. The procedure of condensing 
the unit stiffness matrices before they are combined into the element 
stiffness matrix is valid only if all the subtriangles are coplanar. 
Therefore, if the four corner nodes do not lie on the same plane, i.e., 
if the quadrilateral is warped, errors will occur in the computation of 
element properties. 


If a triangle, rather than a quadrilateral, is used in the program. 



Figure 1.1.2 


a single LCCT^-Q element will be 
used for bending, and a constant 
strain triangle for the membrane 
action. The triangular element 
with its local coordinate system 
is shown in Fig. 1.1.2. 


Triangular Plate Element. 
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The basic element loads consist o£ lateral pressure, thermal loading, 
and gravity loads in the three global coordinate directions. The temper- 
ature is assumed to be uniform through the thickness of the plate, i.e., 
thermal stresses caused by temperature gradients in the direction normal 
to the element are neglected. 

Provision has been made for temperature-dependent material properties. 
The properties of each material used should be listed for a range of 
temperatures covering che entire temperature spectrum of the structure. 
Linear interpolation is used by the computer to calculate the properties 
of each element for the specified nodal temperatures. 

The output from analysis contains the stress resultants referred 
to the coordinates x and y. The positive directions of the stress re- 
sultants are shown in Fig. 1.1.3. For a quadrilateral, the stress re- 
sultants are calculated at node 5 (Fig. 1.1.1) by averaging the values 



Positive Stress Resultants. 



1. 1,4 


for the four sub-triangles. In the case of the triangular element, 
the stress resultants are computed at node 3 in Fig. 1.1.2. 

If the plate elements are used to model a curved surface, the nodal 
rotation about the normal to the surface may be suppressed by the use of 
boundary elements (see Ch. J) , 

The element data deck for each construction code used must begin 

with : 

VI A Element Control Card (515) 


1 

6 

11 

16 

21 25 

MTYPE 

NUME 

NUMMAT 

NUMTC 

KODE j 



JD A D . 


_ ^ 


rvi\ j 




NPAR(l) = Code for element type (MTYPE) . For plate elements use 
NPAR(l) = 6, 

NPAR(2) = Number of elements with the specified construction code (NUME) . 

NPAR(3) = Number of different materials used for the specified construction 
code (NUNIMAT) . 

NPAR(4) = Maximum number of temperatures for which the properties of any 
one material are given (NUMTC) . 

NPAR(5) = Construction code (KODE) . 

The remainder of element data is listed separately for each construction code. 
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I . 2 Construction Code No. 1 

This construction code treats isotropic, homogeneous plate elements. 
The design variable (size of the element) is the plate thickness t. 

The stress ratio redesign formula is obtained from the modified 
Von Mises yield criterion 

where a* and t* are the allowable normal and shear stresses referred to 
the ic and y axis, respectively. Different values of a* may be used for 
tension (a*) and compression (a*) . 

Note that (1.2.1) represents the original , isotropic Von Mises yield 
criterion if a* = a* and x* = a*//3. In that case, the design would be 
independent of the angle 3 in Figs. 1.1.2 and 1.1.3; consequently, 3 will 
only control the printout of the stress resultants. 

The stress ratio formula obtainable from (1.2.1) is the fourth-order 
equation in tVt: 

(t'/t)"^ - a(t'/t)^ ? b(t'/t) - c = 0, (1.2.2) 

where the minus sign preceding b is applicable to the upper surface of 
the plate and the plus sign to the lower surface, and 

a = (N^/N*)^ + (NyN*)2. +• (Ky/N|)^ - Ny^/(N*)^ 

b = 2[(N^/N*)CM./M*) + CNV/N*) (M./M*) + (N../N|) (M^./M|) ] 

- [(N^/N*)(M./M*) + (NVN*)(M^M*)] 

Ay y A 

c = (M^/M*)2 + (M^M*)^ + (MWM^)^ - (M./M*) (M^/M*) . 

A y ^ ^ X 


( 1 . 2 . 3 ) 
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In (1.2.3) we used the notation 


N* = to* 


N* = tx* 
s 


M* = t 


^ a */6 , 


M* = t t */6 


(1.2.4) 


Equation (1.2.2) is solved by the method of successive iterations, 
the iteration equation being 

(r ) - (a + b/r ) (r ) - c = 0 , 

i.e. , 

=[j(a+b/r^) + /j(a+b/r^)^ + c , (1.2.5 

where is the current value of tVt, and r^^^ represents the improved 
value. The starting value of r^ is taken as one, and the number of 
iterations is limited to ten. Note that for pure bending action or pure 
membrane state of stress b = 0, in which case, (1.2.5) becomes an exact 
expression for t’/t. 

The redesign equation is applied only to the point where the stress 
resultants are calculated, namely node 5 for a quadrilateral, and node 3 
for a triangle. No local instability criteria are used with this con- 
struction code. 

VI B Material Properties Data 

A separate data deck, described below, is required for each material 
used for this construction code. The decks do not have to be in the 
sequence of material numbers. 

Material Control Card (215, FIO.O) 


1 

6 

11 20 

N 

NTC 

WT 
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^ = Material number. 

NTC = Number of temperatures for which the properties of this material 
are given. If blank, computer sets NTC = 1. 

WT = Weight-density of the material. 

Material Properties Cards (7F10.0) 

One card for each temperature for which the properties of this 

material are given. Cards must be in ascending order of temperatures. 


1 

11 

21 

31 

41 

51 

61 70 

T 

E 

V 

a 


a* 

c 

T* 


PMATC?) 


PMAT(1)= Temperature for which the material properties are given (T) . 
PMAT(2)= Young’s modulus of elasticity (E) . 

Pt'lAT(3)= Poisson’s ratio (v) . 

PMAT(4)= Coefficient of linear thermal expansion (a). 

PMAT(5)= Tensile strength • 

PMAT(6)= Compressive strength • If blank, computer sets PMAT(6)=PMAT(5) 
PMAT(7)= Shear strength (x*) . If blank, computer sets PI4AT(7)=0 . 577*PMAT(5) 
VI C Element Load Multipliers (4F10.0) 

Five cards shown below: 


1 11 21 31 40 


A 

B 

C 

D 

A 

B 

C 

D 

A 

B 

C 

D 

A 

B 

c 

D 

A 

B 

c 

D 


E^^UL 


Card 

1: 

lateral pressure 

Card 

2: 

thermal loading. 

Card 

3: 

x-gravity . 

Card 

4: 

y-gravity . 

Card 

5: 

z -gravity . 
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EI'4UL = Fractions of. basic element loads (lateral pressure, thermal 
loading and gravity loads in the three global coordinate 
directions) which are to be included in element loads A,B,C 
and D* 

Element Load Multipliers simply define the element loads A,B,C and D. 
Various multiples of these element loads can be added to each structural 
load condition by the use of Structural Load Multipliers (see Sec. D.l). 
IV D Element Data (815, 3'FIO.O) 

One card for each element; cards must be arranged in an ascending 
order of element numbers. 


1 

6 11 16 

21 

26 

31 

36 

41 

51 

61 

71 80 

lEL 

1— 1 

L 

IMAT 

INC 

IDV 

PRESS 

REFT 

FRC 

BETA 


IEC4) H 



lEL = Element number. 

IE = Node numbers I,J,K and L of the four corner nodes (see Fig. 1.1.1). 
For a triangular element use L = 0. 

IMAT = Material number of the element. 

INC = Node number increment in automatic generation of element data 
(see note below). If blank, computer sets INC = 1. 

IDV = Design variable number of element. 

PRESS = Lateral pressure (in the positive z-direction) acting on the 
element. 

REFT = Reference temperature (temperature of the stress-free state) 
of the element. 

FRC = The design variable fraction in (B.1.2): A^ = PiDjj^. Note 

that the size of the element Aj^ is the thickness of the plate. 

If blank, computer sets FRC = 1.0. 


= The angle 3 (in degrees) that defines the x and y axes as 

shown in Figs. 1.1.1 and 1.1.2. Note that the stress resultants 
printed out are referred to the x and y axes. 


BETA 
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Automatic element data generation if there exists a series of 

elements lEL = m, m+1, m+2, , which satisfy the following conditions: 

(a) the node numbers form the sequences 

I = i, i+INC, i+2*INC,..., 

J = j, j+INC, j+2*INC,. . 

K = k, k+INC, k+2*INC,. 

L = ^, A+INC, £+2*INC,..., 

(b) the remainder of 'the data is identical for all elements of the 
series, 

then only the last card of the series will be required. The element data 
for the other elements of the series will be generated automatically. 
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J . BOUNDARY ELEMENTS 


A boundary element is essentially a line element with specified 
extensional and/or rotational spring constant. The element can be 
used for the following purposes: 

(a) modelling of elastic supports, 

(b) enforcement of specified nodal displacements or rotations in 
a given direction, 

(c) computation of support reactions. 

The direction of the boundary element can be specified in two ways, 
as shown in Fig. J.l. 



Figure J.l 

Two Options for Specifying Direction n of 
a Boundary Element. 


In method (a) the direction is determined by the structural node N and 
a second node I. The latter may also be a structural node, or a node 
specially created for the purpose. In the latter case, the degrees of 
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freedom of I should be suppressed on the node, card. In method (b) 
the direction of the element is taken as perpendicular to the lines I-J 
and K-L. Again, the points I,J,K and L may be structural nodes, or 
special points listed on the node cards (with suppressed degrees of 
freedom) . This option is particularly useful when boundary elements 
are used to suppress the normal rotations for thin shells. 

The results of the analysis consist of the axial force and/or 
torque in the boundary element. If the boundary element is an extensional 
spring, the axial force is computed from the formula 

P = k6 , (J.l) 


where k is the spring constant and 6 
(see Fig, J.2 ) . 



Positive Displacements and 
Forces Acting on the Structure 


the computed displacement of node N 

Similarly, the torque in a 
rotational spring obtained 
from 

T = ke , (J.2) 

with T being the torque and 
9 the computed rotation. 


Positive directions of the forces and displacements are defined in Fig. J.2. 
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If a nonzero displacement 5* is prescribed for node N, the following 
force is applied to the node prior to analysis: 


P* = k6* . 


(J-3) 


This applied force will result in the specified displacement only if 
the spring constant k of the boundary element is made much larger than 
the corresponding stiffness of the structure. A specified rotation 0* 
is handled in the analogous fashion, 

A regular elastic support (without a prescribed displacement) is 
obtained by setting 6* = 0 and/or 0* = 0, If, in addition, k is made 
sufficiently large, the boundary element will approximate a rigid 
support. The latter option provides a means of enforcing ’^skewed” 
boundary conditions due to, for example, inclined roller supports. 

Caution must be exercised in the use of boundary elements. Round- 
off errors in the input data and computation of element properties have 
the effect of introducing springs in directions other than the desired 
orientation. The spring constants of these undesired "elements" are 
proportional to the roundoff errors. Although some of the computational 
errors are minimized by the use of double precision arithmetic, it 
is still important not to make the spring constants of the boundary 
elements too large and to specify the directions of the elements with 
great precision. 

VI A Element Control Card (215) 

1 6 10 


MTYPE NUME 


-NPAR(2)- 
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NPAR(l) = Code for element type (MTYPE) . For boundary elements use 
NPAR(l) = 7. 

NPAR(2) = Number of boundary elements (NUME) . 

VI B Element Load Multipliers (4F10.0) 

1 n 21 31 40 

A I B I C I D 
-S EMUL(4) *•“ 

EMUL = Fractions of specified nodal displacements and rotations which 
are to be included in element load cases A>B,C and D. 

Element Load ^^ultipliers define the element loads A,B,C and D. Various 

multiples of these element loads can be added to each structural load 

condition by the use of Structural Load Multipliers as explained in 

Sec, D.l, For boundary elements, the load multipliers are useful if 

some structural load cases include specified nodal displacements, while 

others do not . 

VI C Element Data (915, 5X, 3F10.0) 

One card for each boundary element; cards must be in ascending order 
of element numbers. 


1 

6 11 

16 

21 

26 

31 

36 

41 

46 

51 

61 

71 80 

lEL 

N 1 I 

1 J 1 

K 

1 L 

KD 

|kr 

INCI 

blank 

SD 

1 SR 

TRACE 




IE(5) 


5 ^. 



lEL = Element number. 



= Node number of node N (structural node to which the boundary 
element is attached) . 

= Node numbers, I,J,K,L that define the direction of the element. 

If node I only is used as shown in Fig. J.l(a), set J = K = L = 0. 




J.5 


KD = Code for extensional boundary element. 

KD = 0: extensional spring is not used. 

KD = 1 : extensional spring is used. 

KR = Code for rotational boundary element. 

KR = 0: rotational spring is not used. 

KR = 1: rotational spring is used. 

INC = Node increment used in automatic generation of element data 
(see below) . 

SD = Specified displacement of node N. 

SR = Specified rotation of node N (radians). 

TRACE = Spring stiffness for extensional and/or rotation spring. If 
blank, computer sets TRACE = 10.0**10. 

Automatic element generation if there exists a series of elements which 

satisfy the following conditions: 

(a) the element numbers form the sequence 
MEMB = m, m+1, m+2,...; 

(b) the node numbers of successive elements are 

N = n, n+INC, n+2*INC, . . . , 

I = i, i+INC, i+2*INC, ... 

J = j, j+INC, j+2*INC,...,l or, J = K = L = 0 

K = k, k+INC, k+2*INC,..., > for all elements if node I only 

L = Z, £+INC, £+2*INC,..., is used for direction; 

(c) the remainder of the data is identical for all elements of the 
series ; 

then only the last card of the series will be required. The element 
data for other elements of the series will be generated automatically. 
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